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Abstract 

The work to be presented herein was motivated largely by a desire to improve the understanding of 
oscillatory fluid mechanics inside a Stirling engine. To this end, a CFD project was undertaken at 
Cleveland State University with the goal of accurately predicting the fluid dynamics within an engine or 
engine component. Along with the CFD efforts, a code validation project was undertaken at the 
University of Minnesota. The material covered herein consists of four main parts. In section 1, an 
experimental investigation of a small aspect ratio impinging jet is discussed. Included in this discussion is 
a description of the test facilities and instrumentation. A presentation of the collected data is given and 
comments are made. Next, in section 2, a parallel experimental investigation is presented in which the 
same geometry as that of section 1 is used, but the flow conditions are changed from steady unidirectional 
flow to sinusoidally oscillating flow. In section Two, collected data are presented and comments are 
made. In section 3, a comparison is made between the results of sections 1 and 2, namely, sinusoidally 
oscillating flow results are compared to steady, unidirectional flow results from the same geometry. 
Finally, in section 4, a comparison is made between experimentally collected data (the main subject of 
this work) and CFD generated results. Furthermore, in appendix A, an introductory description of the 
primary measurement tool used in the experimental process — the hot wire anemometer — is given for the 
unfamiliar. The anemometer calibration procedure is described in appendix B. A portfolio of data 
reduction and data processing codes is provided in appendix C and lastly, a DVD and a roadmap of its 
contents is provided in an appendix D. 


1.0 Unidirectional Flow Investigations 

1.1 Introduction 

This unidirectional experimental program was undertaken to complement an oscillatory flow 
investigation conducted at the University of Minnesota. The oscillatory investigation is discussed 
thoroughly in section 2. We defer the description of the motivation behind these experiments until the 
introduction of section 2. The work that is discussed in this thesis began (chronologically) with oscillatory 
flow visualization experiments. It was decided that it would be valuable and important to investigate the 
flow under unidirectional conditions in the same geometry as that of the oscillatory experiments. The 
thought was that the unidirectional case would be less complicated to model with a CFD program (a 
moving boundary would be replaced with a steady state boundary condition). Thus, a series of 
unidirectional experiments were carried out to capture the important features of the flow within the test 
section. The purpose of these experiments was to provide a data set for comparison to CFD generated 
velocity fields. Hot-wire anemometry data were taken and flow visualization was conducted as a standard 
for code validation. The flow geometry was simple, such that it could be easily gridded in a CFD 
program. However, the geometry provided separation and transition zones, shear layers and recirculation 
zones. These characteristics made the flow complex and challenging for CFD computation. 

We comment that the order of experiments that produced this report is as follows: experimental flow 
visualization under oscillatory flow conditions was carried out; this was followed by unidirectional flow 
visualization and hot wire anemometry; finally, oscillatory hot wire anemometiy was conducted. We 
present the results out of chronological order for the following reason: the unidirectional results are easier 
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to understand. We sacrifice order for clarity. In the presentation of unidirectional results that follow, the 
reader will become familiarized with the general flow characteristics within the test section. This 
familiarity will make the presentation of oscillatory results easier to understand and more concise. One 
thing to keep in mind while reading the unidirectional presentation is that the operating points of the 
experiment (under unidirectional flow) were chosen to match the oscillatory case. Also, a large portion of 
the experimental set-up for the unidirectional case is similar, but not identical to, the oscillatory flow case, 
hence, exercise caution when reading section 2. 

Below, the unidirectional experimental process will be outlined and the data collection process will be 
presented. Representative results will be discussed and a summary of all collected data will be given. The 
experimental results will be compared to theoretical predictions and trends will be discussed. 

1.2 Nomenclature 


Symbol 

Description 

Units 

a 

A constant 

[1/s] 

D 

Tube diameter 

[m] 

fit ) 

Generalized function of time 

[various units] 

A(to) 

Fourier transform of / (/) 

[various units] 

F( (0 N ) 

Discrete analogue of F(co) 

[various units] 

Fo 

The NxN Fourier matrix 

[dimensionless] 

gif 

A trial function in z 

[m/s] 

i 

Square root of-1 

[dimensionless] 

N 

An integer 

[dimensionless] 

P 

Pressure 

[N/m 2 ] 

r 

Radial test section coordinate 

[m] 

R 

Tube radius 

[m] 

Re 

Reynolds number 

[dimensionless] 

5 

Disc spacing 

[m] 

t 

Time variable 

[s] 

TI 

Turbulence intensity 

[dimensionless] 

u 

Instantaneous velocity 

[m/s] 

U 

Inviscid velocity 

[m/s] 

u cl 

Tube centerline velocity 

[m/s] 

u j 

Instantaneous velocity 

[m/s] 

u 

Ensemble averaged velocity 

[m/s] 

w 

Boundary layer velocity 

[m/s] 

W 

Invicid velocity 

[m/s] 

X 

Axial test section coordinate 

[m] 

X 

Generalized discrete function of time 

[various units] 

z 

Axial analysis variable 

[m] 

£ 

Error quantity 

[dimensionless] 

<KQ 

Similarity transform function 

[dimensionless] 

R 

Dynamic viscosity 

[N s/m 2 ] 

V 

Kinematic viscosity 

[m 2 /s] 

P 

Density 

[kg/m 3 ] 

X 

Shear stress 

[N/m 2 ] 

CO 

Frequency variable 

[1/s] 

(O n 

Discrete Frequency variable 

[1/s] 

c 

Similarity variable 

[dimensionless] 
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1.3 Experimental Set-Up 


The test hardware consists of the following; a large fan, flexible ductwork, two large plywood flow- 
conditioning boxes, a Laminar Flow Element (LFE), a baffle plate, “polarfleece” quieting sheets, a 
bellmouth and a test section (see fig. 1.1). In figure 1.1, atmospheric air enters the fan (not shown) and is 
pumped through flexible tubing into the first flow conditioning box. This flow conditioning box serves to 
“spread the flow out” and provide an even distribution to the LFE. The LFE kills any large scale eddies 
that may have been caused by the fan. Upon leaving the LFE, the air enters the second flow conditioning 
box and immediately impinges on a baffle plate. This baffle plate destroys the air jet exiting the LFE. The 
flow then passes through the first of two felt-like quieting sheets. Both quieting sheets serve to eliminate 
any eddies that may have remained in the flow. Finally, the flow passes through the second flow 
conditioning box and exits through a bellmouth into the test section. In figure 1.1, two dimensions of 
interest, S and D, are identified. The spacing between the two discs, S, is a variable. The delivery tube 
diameter (or “tube” for short) is denoted by D and is a constant and equal to 216 mm (8.5 in.). For future 
reference, the area inside the tube will be called the tube space and the area between the two discs will be 
called the disc space. It is also to be noted that the flow conditioning boxes are parallelepipeds, whereas 
the laminar flow element and the test section are bodies of revolution (and hence axisymmetric). 

Within the test section, the flow travels through the tube and impinges on the “target” disc 
(see figure 1.2). The flow stagnates, turns radially outward, and exits the test section. An approximate 
time-average picture of the flow inside the test section is presented in figure 1 .2. Flow entering the test 
section (i.e., in the tube) looks generally like a developing pipe flow and is shown in figures 1.33 through 
1.36. For reference, the walls of the test section are made of Plexiglas to allow optical access. Figure 1 .3 
is introduced to give the reader a sense of perspective and size. 

1.4 Data Collection 

Quantitative velocity data are collected with a hot-wire anemometer and probe. The anemometer is a 
model 1053A constant-temperature anemometer made by TSI, Inc. The probe is a TSI model 1210 end 
flow hot-wire sensor (see fig. 1.4). The signal from the anemometer is monitored in one of two ways: for 
low-frequency sampling, the signal is sent directly to an IOTech ADC488 analog-to-digital converter; for 
high frequency sampling, the signal, prior to going to the A/D converter, is run through a D-C offsetting 
device and a Stanford Systems Inc. Dual Channel Filter. Both arrangements are shown in figure 1.5. The 
A/D converter is controlled with C code and is interfaced to a PC via an IEEE 488 bus. The sampling 
stations are shown in figure 1.6. At each sampling station (each radial position), two different types of 
data are collected; a collection of low-sampling-frequency, long-time points and a single, high-sampling- 
frequency point. The collection of low frequency points is used to build an average velocity profile (i.e., u 
vs. x at a fixed r) as the probe is traversed away from the wall in the x-direction (co-ordinate directions 
are shown in figure 1.6). The low-frequency points are collected in a random order (in the x co-ordinate) 
such that any potential temporal variations do not appear as trends associated with x-position. The low- 
frequency sampling rate is 20 Hz, with a sample duration of five minutes. This sampling frequency (i.e., 
20 Hz) was chosen after looking at the results of a high-frequency (2000 Hz) sample. It was observed that 
the Fast Fourier Transform (FFT) of the high-frequency sample shows no significant power contribution 
beyond 20 Hz (shown in figs. 1.13 through 1.32). The concepts of the FFT and “power” are explained 
below. 

For each high sampling frequency point, the wire is approximately 1 mm (0.040 in.) off of the target 
impingement plate. The sampling frequency is 2000 Hz. The choice of frequency says something about 
the expected bandwidth of the signal. The sampling frequency is chosen in accordance with the following 
logic. The flow in the test section has a certain velocity associated with it and, hence, a small packet of 
fluid is convected downstream at this “characteristic” velocity. The packet of fluid (henceforth called an 
eddy) has a certain rotational velocity and spatial size associated with it. Thus, as an eddy flows past the 
anemometer probe, a certain velocity fluctuation will be seen (where the velocity fluctuation can be 
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interpreted as a frequency — that is, the change in velocity due to the eddy looks like a little wave 
propagating past the probe). We choose the sampling frequency to capture the smallest eddies of interest 
traveling at the characteristic fluid velocity. We shall see during the discussion of power-spectra plots that 
this choice of sampling frequency is quite reasonable. The duration of the high frequency sample is 
32.768 seconds. The sampling frequency-time combination gives a total of 65536 = 2 16 data points which, 
by virtue of the “power of 2” size of the sample, optimizes the Fourier transform (to be discussed shortly). 

Flow visualization equipment consists of a Nikon FE2 35 mm camera with a Tamaron SP 1 :2.5 by 
90 mm lens. A schematic of the set up is shown in figure 1.7. In all cases, there is a red filter screwed on 
the lens; the filter serves to cut all non-red light from the photographs. The camera is loaded with Tri-X 
ASA 125 black and white film and mounted on a tripod with the camera shutter release set to “bulb.” 

Bulb is a setting that holds the camera shutter open as long as the shutter trigger is depressed. The camera 
is aligned normally with the x-axis of the test section and is focused such that the x-axis is in the focal 
plane. The focal plane is illuminated with a LASER sheet. The LASER source is a Melles-Griot Fle-Ne 
LASER that emits red light at a frequency of 632.8 nm. A glass rod is placed in the path of the beam to 
create the sheet. Particles are introduced into the flow with a Sage Action Inc. Model 33 Multi Head 
Bubble Generator. This device combines compressed air, helium and a soap solution to make neutrally- 
buoyant, variable-sized bubbles. The bubble size is in general about 0.8 mm (1/32 in.). The bubbles, 
which are injected just downstream of the second quieting screen (see fig. 1.1), travel with the flow into 
the test section. When the bubbles arrive in the disc space, some of them are illuminated in the laser sheet. 
The product is a real-time, two-dimensional view of the flow. To capture a picture of the flow, the camera 
is set-up as described above and the shutter is opened for five minutes. This method of photography 
produces a time average picture (shown in figs. 1 .8 through 1 . 1 1 ) of the flow. As a note, due to the long 
duration of the photo, all of this work is carried out in a blackened room. 

The cases of interest are as follows: 

1. Re = 7600, S= 127 mm (5 in.) 

2. Re = 7600, S = 54 mm (2.125 in.) 

3. Re = 17700,5= 127 mm (5 in.) 

4. Re = 17700, S= 54 mm (2.125 in.) 

where the Reynolds number is based on tube diameter, D, and tube core exit velocity. 

1.5 Results 


1.5.1 Flow Visualization 

The flow visualization scheme produced several photographs. Scanned versions are shown in 
figures 1 .8 through 1.11. The photographs give a rough idea of where transition to turbulence along the 
target wall takes place. To locate transition zones, one looks at the flow adjacent to the target plate. While 
moving the eye radially outward along the near-wall flow, one observes that the streaks start to wiggle. 
The presence of “wigglyness” implies that the flow has begun to transition. The location of the onset of 
transition, as revealed by the flow visualization, is pointed out with white arrows in figures 1.8 through 
1.11 and is summarized below. As a note, in figures 1 .8 through 1.11, two sets of arrows appear. The 
green (or gray if printed in black and white) arrows will be discussed below. 

Re = 17700, S = 127 mm (5 in.) 

Figure 1.8 

Normalized Radial Transition Location ( r/D ) = 1 

Re = 7600, S = 127 mm (5 in.) 

Figure 1.9 

Normalized Radial Transition Location (r/D) ~ 1 
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Re = 17700, S = 54 mm (2.125 in.) 

Figure 1.10 

Normalized Radial Transition Location ( r/D ) = 0.8 

Re = 7600, S = 54 mm (2.125 in.) 

Figure 1.11 

Normalized Radial Transition Location (r/D) ~ 1 

One thing that the reader may note about the flow visualization pictures it that in some cases (figs. 1.8 
and 1.10) the “top” of the photo is more heavily laden with bubbles, whereas in figures 1 .9 and 1.11 the 
bottom is more bubble-dense. This is because the bubbles are not exactly neutrally buoyant. The photos 
were taken on different days, and, as such, the bubble machine was turned off and on. It is difficult to 
exactly control the buoyancy of the bubbles; hence, sometimes they are slightly heavy and fall towards 
the floor in the plenum upstream of the test section where the velocities are very low. On other days they 
rise. The net result is the phenomenon that appears in the photographs. 

1.5.2 Fourier Analysis and Spectral Techniques 

As was mentioned in the data collection section, a high-sampling rate point was taken at each radial 
location for all four cases. We use these signals to help determine the location of transition to turbulence 
in the boundary layer. 

To begin, the Fourier transform of a continuous signal, and its implications, will be described. The 
discrete case, which is inherent to data collection, is discussed thereafter. 

The Fourier transform of a continuous signal can be defined as: 


F( co)= \f(t)e‘ m dt 


where f(t ) is a real-time signal (i.e., a continuous and bounded function). Integration against a complex 
exponential takes the time signal to “frequency space.” The idea is that at a fixed frequency, (0, the 
product of the time signal and the complex exponential will “pick out” the oscillations of the time signal 
that correspond to the given frequency, to. In short, we now have a function that describes the frequency 
content of our original time signal. 

The description of how this “picking out” process works is mathematically deep and the curious 
reader is directed to Pinsky 2 (or your favorite text on Fourier analysis). This integration process produces 
a (usually) continuous function in “frequency space.” The new function is (almost always) complex, and, 
as such, requires some kind of processing. To this end, the function F(co)is pointwise squared (i.e., we 

9 

take the product F(co)xF(co)at each co). We define this new function F~(c o) as “power” or “energy.” As 
a note, F(co) carries the units (in the present discussion) of length/frequency. Thus, the square of this 
combination of units produces a quantity that is similar in form to the idea of kinetic energy, hence the 
name “energy” or “power.” 

Data collection is inherently discrete, hence the continuous signal f(t) is sampled at some frequency 
co. Suppose N points are collected. Then the data set can be represented as a column vector called X, 
where X has a length of N. 
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where the V|/'s are defined as the N roots of \|/' V =1. The Discrete Fourier Transform (DFT) and the 
power associated with it are defined now: 


F(c 0 N ) = FoX 

and 

F~(co n ) = F(q) n )F h (co N ) 

We note that there are optimization algorithms associated with the above definition of the DFT. We 
do not care to elaborate on the details but rather refer the reader to the references (Flayes, 1 Strang 4 ). We 
note further that these algorithms take the DFT to the Fast Fourier Transform (FFT). The FFT depends on 
samples (vectors) of length 2 J (recall that our sample length, N, was chosen to meet this requirement). 

Figures 1.13 through 1.32 present two plots each; one is the actual time signal, X, and the second is 

— 9 

the “power,” F ( co N ) . We use spectral (i.e., high-frequency Fourier transformed) data to determine the 
location of transition onset. The general method is as follows: one observes the frequency domain 
velocity signals (i.e., power) at the radial sampling positions r/D = 0.1 18, 0.471, 0.588, 0.824, and 0.941. 
One can see that as r/D increases, the power signal first flattens, then develops a “hump” in the 
neighborhood of 2 ITz < co < 20 FIz (figs. 1.13 through 1.32). The range of frequencies associated with 
transition is shown in figure 1.12. We associate transition with a power rise in the vicinity of these 
frequencies. The logic is that a fluctuation power increase over a certain frequency range corresponds to 
(in the physical situation) the emergence of coherent waves being carried past the probe. The presence of 
coherent waves implies the onset of turbulence. 

This process was carried out for all four cases of different Reynolds numbers and disc spacings. It 
was found that there was general agreement with the flow visualization results, with one exception; the 
transition points were generally at a slightly smaller upstream radial position than the transition onset 
locations that the flow visualization data indicated. In other words, the hot-wire revealed the development 
of turbulence to a higher level of sensitivity than that of the flow visualization. The insight that the 
spectral technique provides is shown with green arrows (or gray if printed in black and white) in contrast 
to the flow visualization in figures 1 .8 through 1.11 and is given quantitatively below. 
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Re = 17700, S = 54 mm (2.125 in.) 

Figures 1.13 to 1.17 and 1.10 

Normalized Radial Pre-Transition Location (r/D) = 0.588 
Normalized Radial Post-Transition Location ( r/D ) = 0.824 

(On fig. 1.10, transition is identified as r/D = 0.7 for the sake of comparison to flow visualization 
data, we can only say that it is between 0.588 and 0.824.) 

Re = 17700, S = 127 mm (5 in.) 

Figures 1.18 to 1.22 and 1.8 

Normalized Radial Transition Location (r/D) = 0.941 

Re = 7600, S = 54 mm (2.125 in.) 

Figures 1.23 to 1.27 and 1.11 

Normalized Radial Transition Location (r/D) = 0.941 

Re = 7600, S= 127 mm (5 in.) 

Figures 1.28 to 1.32 and 1.9 

Normalized Radial Transition Location (r/D) = 0.824 


1.5.3 Ensemble Averaged Velocity and Turbulence Intensity Profiles 

Velocity and turbulence intensity profiles were generated using long-time, low-frequency sampling. 
In the disc space, these data offer the primaiy test of code validation. To seed the computations, tube 
velocity and turbulence intensity profiles are used as a boundaiy condition for CFD calculations. 
Turbulence intensity is defined in this case as: 


j t( u j~ u f 

77 = —1 — 

u 1 N - 1 

Where N is the number of sampled points, u is the local velocity averaged over all of the sampled data 
points and u j is the instantaneous velocity of the /th point. 

The boundary data are (i.e., tube velocity and 7/ profiles) shown in figures 1.33 through 1.40. The 
quantity being measured is u , the average velocity in the axial direction (the assumption is that in the 
tube space, the main component if the velocity is in the axial direction). The axial location of the probe in 
the tube is (as is shown in fig. 1.6) 127 mm downstream of the test section entrance. There is apparently a 
slight skewing of the flow which approaches the bellmouth entrance to the test section. The reader might 
notice that the profile shapes are not exactly symmetric (in terms of velocity and 77). One other thing to 
note is that there are considerably more data points taken in the radial range 0 < r/D <0.5 than over the 
range - 0.5 < r/D < 0. Emphasis was given to one side because the team at CSU used only one half of the 
measured domain for boundary conditions. (They made the assumption that the flow is axisymmetric, 
hence it is two-dimensional. Thus, assuming that the flow is two-dimensional, they reduced the 
computation domain to cover one -half of the tube). The sparse data on the other side are given to 
document symmetry. 

Turbulence intensity profiles in the disc space are presented in figures 1 .41 through 1 .44. These 
turbulence intensity profiles are useful in the sense that they indicate where the hot-wire is collecting 
“reasonable” data. In other words, in figures 1.41 through 1.44, regions exist where 77 values exceed 
25%. Measured values of velocity rise above actual values and measured values of 77 drop below actual 
values as 77 rises above 25%. 
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Average velocity profiles are presented in figures 1 .45 through 1 .48. Using the 77 data, one can divide 
the velocity profiles into two regions; a region where average velocities may be considered valid for code 
validation purposes and a second region where the profiles are not accurate. In figures 1.45 through 1.48, 
entire average velocity profiles are presented. In regions of large 7/ the points along the profile are 
colored black. 

Using the data from the average velocity and 7/ profiles, we are able to distinguish between the near- 
wall impinging-jet-type flow and the reverse flow far off the target plate and identify the shear layer that 
separates the two. We defer further discussion of this topic to section 1.7. 

1.6 Average Velocity profiles and a Description of the Shear Layer 

Of interest was the documentation of the shear layer (or recirculation zone) shape (see fig. 1.52). To 
accomplish this task, the peak velocity at each normalized radial location was observed. This peak 
velocity corresponds with a particular normalized distance from the target wall. For a clarification of this 
process, we refer the reader to figure 1.49. The data for figure 1.49 came from figure 1.48. In figures 1.45 
through 1 .48 ensemble averaged velocity profiles as a function of normalized distance (x!S) from the 
target wall are presented. Figures 1.45 through 1.48 present five profiles, where the profiles correspond 
to the normalized radial positions (r/D) 0.1 18, 0.471, 0.588, 0.824, and 0.941. At each operating point (as 
functions on RE and S), the velocity profiles corresponding to r/7) = 0.588, 0.824, and 0.941 were 
inspected. From each of these cases, (12 in all) the maximum value of velocity was chosen. This 
collection of maximum velocities was then plotted and is shown in figure 1.50. 

Figure 1.50 reveals that the shape of the shear layer seems to be independent of Reynolds number. 
However, one can see that curves of equal disc spacing seem to follow each other. This implies that the 
shape of the shear layer with radial location, r/D, depends primarily on dimensionless disc spacing, S/D. 

1.7 Experimentally Collected Data vs. Theoretical Results 

1.7.1 An “Exact” Solution of the Navier-Stokes Equations 

As a matter of curiosity, the time-averaged velocity profiles collected via low sampling rate 
anemometry are compared to an “exact” solution of the Navier-Stokes equations. The method of 
comparison is given below. In the viscous region, a theoretical boundary layer analysis is compared to 
experimentally measured shear stresses. 

Before we begin, we make a comment; we note that in figures 1 .33 through 1 .36, the delivery tube 
centerline velocities, u cl , are, respectively, 1.23, 1.15, 0.50, and 0.52 m/s. With the definition: 


Re=B^ 

and air at Standard Temperature and Pressure (STP), we find that these centerline velocities correspond to 
Reynolds numbers of, respectively, 16422, 15354, 6676, and 6943. Thus, the previously referenced 
Reynolds numbers of 17700 and 7600 were only approximate. This difference is ineffectual in the sense 
that code validation was based on the velocity profiles which yielded the correct Reynolds numbers rather 
than the nominal values. In the analysis that follows, the value of the Reynolds number at the tube 
centerline is important, so we make the following formal statement: the cases of interest are changed to: 

Case 1. Re = 16422, S = 54 mm (2.125 in.) (corresponding to fig. 1.33) 

Case 2. Re = 15354, S = 127 mm (5 in.) (corresponding to fig. 1.34) 

Case 3. Re = 6676, 5 = 54 mm (2.125 in.) (corresponding to fig. 1.35) 

Case 4. Re = 6943, S = 127 mm (5 in.) (corresponding to fig. 1 .36) 
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The laminar solution of the Navier-Stokes equations in an axisymmetric stagnation flow is well 
known and is presented in sources such as “Boundary Layer Theory” by Schlichting. 3 In the inviscid 
region (see fig. 1.51), potential flow theory yields the following equations for the U and W velocities: 


U = ar 

(1) 

W = -2 az 

(2) 


where a is a constant and is related to the “strength” of the jet. Close to the impingement target surface 
there will be viscous effects. In that region, the following form of the Navier-Stokes equations govern the 
fluid dynamics: 


du du 1 dp 

u — + w — = — +V| 

or dz p dr 


f d 2 u 1 du 


dr 


- + 


dr 


- + - 


d 2 u) 


ar 


dw dw 1 dp 

u — +w— = — + V| 

dr dz p dz 

du u dw 

37 + 7 + aP° 


( d 2 w 


i du a w 
+-— +■ 


dr 2 r dr dz" 


With the boundary conditions: 


u(r, 0) = w(r,0) = 0 
u(r,z -> oo) = u 


Introducing the transformations (with primes denoting differentiation with respect to z): 

u = rg'(z) 
w = -2 g(z) 


into (3a) and (3b) yields: 


g' 2 -2gg'=a 2 +vg"' 


and 


g(0) = g'(0) = 0 
g'(z^°°) = a 


Now a similarity transformation is introduced: 



g(z) = (KOa/ov 

Putting this transformation into (4a) and (4b) gives: 


(3a) 


(3b) 


(4a) 


(4b) 


(5a) 

(5b) 
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(6a) 


4>' ' ' +2cf>())' ' — cf)' 2 +1 = 0 

4>(0) = cj)' (0) = 0 (6b) 

<t>'(C-»~) = l (6c) 


(Here, primes denote differentiation with respect to £.) 

The solution to equation (6a) subject to ((6b) and (6c)) is presented in Schlichting, 3 pp. 90-91. 

We now find an exact solution for wall shear stress. Regarding equation (5a), we form the differential 
and write: 


d_ = ± Jv 
dCdzia 

From the algebra required to go from equations (4) to equations (6), we know that: 


f= 


u 

U 


Combining equations (1), (7) and (8) and introducing viscosity yields: 


T(Z) = 



( 7 ) 


( 8 ) 


If we are interested in the shear stress at the wall, we have: 


x 


W 



( 9 ) 


We now pause to review. So far we have developed a relationship for wall shear stress as a function 
of radial position along the target plate. We also state that this solution is exactly valid for the geometry of 
figure 1.51. 

The method of comparison is as follows; assume that in the neighborhood of the stagnation point (see 
fig. 1.52) the exact solution closely models the experimental set-up. Thus, we can use (9) to predict wall 
shear stress close to the tube centerline (we suspect that away from the centerline, the flow, and 
specifically the boundary layer, will feel the influence of the shear layer and recirculation bubble above 
it). The next assumption is about the inviscid flow region. We assume that the tube exit velocity at the 
centerline u cl , can be specified at a distance 5 from the target plate. With these thoughts in mind, we turn 
to equation (2): 


W = -2 az 


Inserting our assumptions we have: 


u c l = -2 aS =>l a I 


«d _ Re P v 
25 2SD 


( 2 ) 


( 10 ) 
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Where the Reynolds number is based on tube diameter, D and the tube centerline exit velocity u c [ . 
Putting (9) and the result of (10) together gives: 


[iV(0) 

Re» 

<N 

Cl 

00 

i 

S 3 


( 11 ) 


where R is the tube radius. The known numerical values are: 

R = 0.10759 m (4.25 in.) 

p= 1.1614 kg/m 3 
jj. = 184.6 xlO -7 N • 5 / m 2 
<|>"(0) = 1.3 120 [dimensionless] 


This gives (as an approximation): 


1.364 xlO -9 . 

Re!> 

\ 

S 3 


( 12 ) 


Near the stagnation point we expect that the above equation will be in close agreement with measured 
shear stress results. We hope that reasonable agreement will extend from the stagnation point to a radial 
location equal to the radius of the tube. 

We compare shear stress data at radial locations r = 102 mm (4 in. or r/D = 0.473) and r = 127 mm 
(5 in. or r/D = 0.588). We presume that in this region, and for smaller values of r, the effect of the 
difference in geometry between that of figure 1.51 and the actual geometry of figure 1 .52 will be minimal. 
We note that the data-collection process was not originally intended for shear stress measurements in that 
we did not focus on near- wall measurements. As such, some of the data give a reasonable approximation 
of the gradient (and consequently the wall shear stress) whereas, others do not. One can get a sense of the 
accuracy of the shear stress by looking at the near wall velocity gradients in figures 1 .45 through 1 .48 . 
Only the “good” data will be presented. Finally, we note that these data represent only the lower bounds 
to the actual values of wall shear stress. (The closest data point was 1 mm from the target wall, thus, the 
gradient at 1 mm was assumed to extend all the way to the wall.) 

Regarding figure 1 .53, one can see that the experimental values are in agreement with the theoretical 
values and the trends with S and Re D (actually u cl ) are captured reasonably well. 

1.8 Conclusions 

The flow field in a semi-constrained, steady, impinging jet was investigated. The flow was 
documented with flow visualization and hot wire anemometiy. The flow visualization results served to 
provide a qualitative characterization of the flow field. From it, the general shapes of the shear layer and 
wall jet could be identified. Moreover, from the flow visualization results, approximate transition points 
in the wall jet could be identified. With the knowledge gained from the flow visualization, we turned to 
hot wire anemometry coupled with Fourier theory to identify wall jet transition zones. Low frequency hot 
wire samples allowed for the generation of average velocity and 77 plots that are useful for code 
validation purposes. From the low frequency hot wire data, plots of the shear layer shape were made. It 
was found that the shear layer is more or less independent of Reynolds number and strongly dependant on 
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disc spacing. Finally, approximate values of wall shear stress were computed. These results, by 
themselves, are useful for code validation; however, it was found that when compared to analytical 
predictions of shear stress, the results were very favorable; a much more general and useful result. 


1.9 Error Analysis 

A hot wire calibration error analysis was conducted and is described in detail in appendix B. It 
revealed that the hot wire anemometer has a variable uncertainty associated with it. The uncertainty in 
hotwire measurements is presented in appendix B, figure B.5. This figure shows that as fluid velocities 
drop below 20 cm/s there is a marked increase in error. We note, however, that even at low velocities (sub 
20 cm/s) the measurement is very repeatable; thus, we do not expect a large stochastic error. For a 
detailed explanation the curious reader is referred to appendix B. 

Outside of anemometer errors, there are errors associated with the positioning of the hot wire probe. 
The probe saw the flow in two different regions of the test section; the “tube” and the “disc space.” 

Within the tube, the probe was allowed to move only radially. However, in the disc space, the probe was 
moved both radially and axially. Each of these positioning schemes has its associated errors, which we 
now state. 

Within the tube, for any radial position, there is a normalized error (error/tube diameter) of: 

£ tube = ±0.00094 

In other words, for any radial position within the tube (e.g., figs. 1 .33 through 1 .40), the plotted or stated 
values may have the above deviation in radial position. 

Within the disc space, for any axial location (i.e., normalized distance from the wall) there is 
normalized error (error/disc spacing) of: 

£ disc,axi, 2 . 125 - = ±0-0047 (for the 54 mm (2. 125 in.) case) 

£ disc,axi,5 " = ±0-002 (for the 127 mm (5 in.) case) 

In other words, for any axial position within the disc space (e.g., fig. 1.47), the plotted or stated values 
may have the above deviation in axial position. 

Finally, the probe was moved radially within the disc space. The normalized error (error/tube 
diameter) associated with radial position within the disc space is stated as: 

^ disc, radial ~ ± 0-012 

This type of error may appear in statements such as rID = 0.941 . 

For completeness, all of the above error statements are based on a 95% confidence interval. 
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1.11 Figures 



Figure 1.1. — A schematic of the experimental set-up. Flow enters on the right and passes through 
the various flow conditioning elements. As a note, the flow conditioning boxes 
are parallelepipeds, whereas the laminar flow element and the test section are bodies of 
revolution (and hence axisymmetric). 



Figure 1. 2. — An approximate picture of the flow within the test section. 
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isc major diameter = 1 ,22m 




12 7 mm ( 50) Hoi Film 
9.5 mm (.38) Hoi W» r e 

38 mm (150) — 

I 

— I 


L 3.2 mm (125) Dia. 



i 


4.6 mm (.tB) Dia. 


Figure 1.4. — A TSI model 1210 hot-wire probe. 
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Low-Frequency Sampling High-Frequency Sampling 

Figure 1.5. — A schematic of the data collection hardware. 



r = 25.4mm ( 1 ") 


r = 102mm (4") 
r = 127 mm ( 5 ") 


Axis of Symmetry 


r = 1 78mm (7") 
r = 203mm (8”) 


Tube Sampling Location 


Figure 1.6. — A cross-section of the test section. Sampling stations are indicated with 
arrows. The “tube traverse location” was 127 mm (5 in.) downstream of the test 
section entrance. The tube traverse was in the “r” direction. All of the traverses in 
the disc space were in the “jc” direction. This is also the view that the camera sees in 
flow visualization experiments. Essentially, the x-axis is in focus (i.e., r = 0). 
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Flow Visualization Set-Up 


(side view) 

Figure 1.7. — Flow visualization hardware set up schematic. 
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Figure 1.8. — Re = 17700, S = 127 mm (5 in.). Flow visualization indicates 
that transition occurs at approximately r/D = 1 (labeled with white 


arrows). Spectral data indicate that transition occurs at approximately r/D 
= 0.941 (labeled with green arrows, or gray if printed in black and white). 
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Figure 1 .9. — Re = 7600, S = 127 mm (5 in.). Flow visualization indicates that 
transition occurs at approximately r/D = 1 (labeled with white arrows). 
Spectral data indicate that transition occurs at approximately r/D = 0.824 
(labeled with green arrows, or gray if printed in black and white). 
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Figure 1.10. — Re = 17700, S = 54 mm (2.125 in.). Flow visualization 
indicates that transition occurs at approximately rID = 0.8 (labeled with 
white arrows). Spectral data indicate that transition occurs at approximately 
r/D = 0.706 (labeled with green arrows, or gray if printed in black and 
white). 
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Figure 1.11 . — Re = 7600, 5 = 54 mm (2.125 in.). Flow visualization indicates 
that transition occurs at approximately r/D = 1 (labeled with white arrows). 
Spectral data indicate that transition occurs at approximately r/D = 0.941 
(labeled with green arrows, or gray if printed in black and white). 
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time domain signal 



Frequency, Hz 


Figure 1.12. — A typical FFT of a velocity-time signal. The “power” and the 
corresponding frequencies in the oval are associated with transition to turbulence. 




Frequency, Hz 


Figure 1.13 . — A picture of the velocity-time and velocity-FFT signals for the 
case RE = 17700, S = 54 mm (2.125 in.), r/D = 0.118. The distance from the 
impingement target plate = 1 mm (0.040 in.). Sampling frequency = 2000 FIz. 
Cutoff frequency = 750 FIz. 
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RE =17700, S = 2.125", r/D = .471 




Frequency, Hz 

Figure 1 .14. — A picture of the velocity-time and velocity-FFT signals for the 
case RE = 17700, S = 54 mm (2.125 in.), r/D = 0.471. The distance from the 
impingement target plate = 1 mm (0.040 in.). Sampling frequency = 2000 FIz. 
Cutoff frequency = 750 Hz. 


RE = 17700, S = 2.125", r/D = .588 




Frequency, Hz 

Figure 1.15 . — A picture of the velocity-time and velocity-FFT signals for the 
case RE = 17700, S = 54 mm (2.125 in.), r/D = 0.588. The distance from the 
impingement target plate =1 mm (0.040 in.). Sampling frequency = 

2000 Hz. Cutoff frequency = 750 Hz. 
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RE = 17700, S = 2.125", r/D = .824 


2 - 



0 5 10 15 20 25 30 35 

Time, s 



Frequency, Hz 

Figure 1.16. — A picture of the velocity-time and velocity-FFT signals for the 
case RE = 17700, S = 54 mm (2.125 in.), r/D = 0.824. The distance from the 
impingement target plate =1 mm (0.040 in.). Sampling frequency = 2000 FIz. 
Cutoff frequency = 750 FIz. 


RE = 17700, S - 2.125", r/D = .941 

1 1 



0 5 10 15 20 25 30 35 

Time, s 



Frequency, Hz 

Figure 1.17. — A picture of the velocity-time and velocity-FFT signals for the 
case RE = 17700, 5 = 54 mm (2.125 in.), r/D = 0.941. The distance from the 
impingement target plate = 1 mm (0.040 in.). Sampling frequency = 

2000 Hz. Cutoff frequency = 750 Hz. 
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10 10 10 
Frequency, Hz 


Figure 1.18. — A picture of the velocity-time and velocity-FFT signals for the 
case RE = 17700, 5 = 127 mm (5 in.), r/D = 0.1 18. The distance from the 
impingement target plate = 1 mm (0.040 in.). Sampling frequency = 

2000 Hz. Cutoff frequency = 750 Hz. 


RE = 17700, S = 5", r/D = .471 




Frequency, Hz 


Figure 1.19. — A picture of the velocity-time and velocity-FFT signals for the 
case RE = 17700, 5=127 mm (5 in.), r/D = 0.47 1 . The distance from the 
impingement target plate = 1 mm (0.040 in.). Sampling frequency = 

2000 Hz. Cutoff frequency = 750 Hz. 
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RE = 17700, S = 5", r/D = .588 




- 2-101234 

10 10 10 10 10 10 10 


Frequency, Hz 


Figure 1.20. — A picture of the velocity-time and velocity-FFT signals for the 
case RE = 17700, S = 127 mm (5 in.), r/D = 0.588. The distance from the 
impingement target plate = 1 mm (0.040 in.). Sampling frequency = 

2000 FIz. Cutoff frequency = 750 Hz. 


RE = 17700, S = 5", r/D = .824 



0 5 10 15 20 25 30 35 

Time, s 



Frequency, Hz 


Figure 1.21 . — A picture of the velocity-time and velocity-FFT signals for the 
case RE = 17700, S = 127 mm (5 in.), r/D = 0.824. The distance from the 
impingement target plate = 1 mm (0.040 in.). Sampling frequency = 

2000 Hz. Cutoff frequency = 750 Hz. 
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RE - 17700, S . 5", r/D - .941 




Figure 1.22. — A picture of the velocity-time and velocity-FFT signals for the 
case RE = 17700, S = 127 mm (5 in.), r/D = 0.941 . The distance from the 
impingement target plate = 1 mm (0.040 in.). Sampling frequency = 

2000 Hz. Cutoff frequency = 750 Hz. 


Re = 7600, S = 2.125", r/D = .118 
Distance from impingement plate = .040" 




- 2-10 1 2 3 4 
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Frequency, Hz 


Figure 1.23. — A picture of the velocity-time and velocity-FFT signals for the 
case RE = 7600, S = 54 mm (2.125 in.), r/D = 0.1 18. The distance from the 
impingement target plate = 1 mm (0.040 in.). Sampling frequency = 

2000 Hz. Cutoff frequency = 750 Hz. 
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Re = 7600, S = 2.125", r/D = .471 
Distance from impingement plate = .040" 




Figure 1.24. — A picture of the velocity-time and velocity-FFT signals for the 
case RE = 7600, S = 54 mm (2.125 in.), r!D = 0.47 1 . The distance from the 
impingement target plate = 1 mm (0.040 in.). Sampling frequency = 

2000 Hz. Cutoff frequency = 750 Hz. 


Re . 7600, S . 2.125", r/D . .588 
Distance from impingement plate = .040“ 




Frequency, Hz 


Figure 1.25. — A picture of the velocity-time and velocity-FFT signals for the 
case RE = 7600, S = 54 mm (2.125 in.), r/D = 0.588. The distance from the 
impingement target plate = 1 mm (0.040 in.). Sampling frequency = 

2000 Hz. Cutoff frequency = 750 Hz. 
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Re = 7600, S = 2.125", r/D = .824 
Distance from impingement plate = .040" 




Figure 1.26. — A picture of the velocity-time and velocity-FFT signals for the 
case RE = 7600, S = 54 mm (2.125 in.), r/D = 0.824. The distance from the 
impingement target plate = 1 mm (0.040 in.). Sampling frequency = 

2000 FIz. Cutoff frequency = 750 Hz. 

Re . 7600, S - 2.125". r/D . .824 
Distance from impingement plate = .040" 




Frequency, Hz 


Figure 1.27. — A picture of the velocity-time and velocity-FFT signals for the 
case RE = 7600, S = 54 mm (2.125 in.), r/D = 0.941. The distance from the 
impingement target plate = 1 mm (0.040 in.). Sampling frequency = 

2000 Hz. Cutoff frequency = 750 Hz. 
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Frequency, Hz 


Figure 1.28. — A picture of the velocity-time and velocity-FFT signals for the 
case RE = 7600, S = 127 mm (5 in.), r/D = 0.1 18. The distance from the 
impingement target plate = 1 mm (0.040 in.). Sampling frequency = 

2000 Hz. Cutoff frequency = 750 Hz. 




Frequency, Hz 


Figure 1.29. — A picture of the velocity-time and velocity-FFT signals for the 
case RE = 7600, 5=127 mm (5 in.), r/D = 0.47 1 . The distance from the 
impingement target plate = 1 mm (0.040 in.). Sampling frequency = 

2000 Hz. Cutoff frequency = 750 Hz. 
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Figure 1 .30. — A picture of the velocity-time and velocity-FFT signals for the 
case RE = 7600, S = 127 mm (5 in.), r/D = 0.588. The distance from the 
impingement target plate = 1 mm (0.040 in.). Sampling frequency = 

2000 FIz. Cutoff frequency = 750 Hz. 


RE = 7600, S = 5", r/D = .824 
0.25 1 1 1 1 
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Figure 1.31 . — A picture of the velocity-time and velocity-FFT signals for the 
case RE = 7600, S = 127 mm (5 in.), r/D = 0.824. The distance from the 
impingement target plate = 1 mm (0.040 in.). Sampling frequency = 

2000 Hz. Cutoff frequency = 750 Hz. 
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Figure 1 .32. — A picture of the velocity-time and velocity-FFT signals for the 
case RE = 7600, S = 127 mm (5 in.), r/D = 0.941 . The distance from the 
impingement target plate = 1 mm (0.040 in.). Sampling frequency = 

2000 Hz. Cutoff frequency = 750 Hz. 
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Figure 1.33. — Average velocity profde in the tube space for the case Re = 17700, 

S = 54 mm (2.125 in.). 
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; 1.34. — Average velocity profile in the tube space for the case Re = 17700, 
127 mm (5 in.). 
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Figure 1.35. — Average velocity profile in the tube space for the case Re = 7600, 

S= 54 mm (2.125 in.). 
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; 1.36. — Average velocity profile in the tube space for the case Re = 7600, 
127 mm (5 in.). 
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Figure 1.37. — Turbulence intensity in the tube space for the case Re = 17700, 

S = 54 mm (2.125 in.). 
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: 1.38. — Turbulence intensity in the tube space for the case Re = 17700, 
127 mm (5 in.). 
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Figure 1.39. — Turbulence intensity in the tube space for the case Re = 7600, 
S = 54 mm (2.125 in.). 
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Figure 1.40. — Turbulence intensity in the tube space for the case Re = 7600, 
S = 127 min (5 in.). 
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Figure 1.41 . — Turbulence intensities in the disc space for the case Re = 17700, 
S = 54 mm (2.125 in.). 
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Figure 1.42. — Turbulence intensities in the disc space for the case Re = 17700, 
S= 127 mm (5 in.). 
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Figure 1.43. — Turbulence intensities in the disc space for the case Re = 7600, 
S = 54 mm (2.125 in.). 
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Figure 1.44. — Turbulence intensities in the disc space for the case Re = 7600, 
S = 127 mm (5 in.). 
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Figure 1.45. — Average velocities in the disc space for the case Re = 17700, S = 54 mm 
(2.125 in.). Portions of the velocity profiles in regions of high turbulence intensity are 
colored black. 
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Figure 1.46. — Average velocity in the disc space for the case Re = 17700, S = 127 mm 
(5 in.). Portions of the velocity profiles in regions of high turbulence intensity are 
colored black. 
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Figure 1.47. — Average velocity in the disc space for the case Re = 7600, S = 54 mm 
(2.125 in.). Turbulence intensity never exceeds 25% in these profiles. 
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Figure 1.48. — Average velocity in the disc space for the case Re = 7600, S = 127 mm 
(5 in.). Portions of the velocity profiles in regions of high turbulence intensity are 
colored black. 
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Figure 1.49. — The figure above illustrates the method of determining wall 
jet thickness. 
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Figure 1.50. — The figure above shows the development of the shear layer for the 
Reynolds number - disc spacing cases of interest. 



Figure 1.51 . — This is the geometry that is applicable to the 
analysis covered in Schlichting. The impingement plate should 
extend to infinity. 
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Figure 1.52. — The test section geometry at the University of 
Minnesota. We expect that equation (12) applies in the 
neighborhood of the stagnation point. The recirculation zones 
have been observed via flow visualization and are included to 
give the reader a better understanding of the flow field. 
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Figure 1.53. — Experimentally gathered shear stress values vs. theoretical predictions. 
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2. Oscillatory Flow Investigations 

2.1 Introduction 

This study addresses the need for validation of multi-dimensional, oscillatory flow techniques that 
can be employed in codes used for Stirling engine design. There is concern about the accuracy of 
modeling certain regions, such as the compression space, of a Stirling engine. This is addressed by 
isolating an important region of the engine (e.g., the compression space) and identifying the important 
features of the flow there. Once an area of interest is chosen, experiments are conducted to characterize 
the flow, compute the flow in those experiments, identify weaknesses in the computation and modify the 
models to improve their ability to simulate the experiments. The result is that Stirling engine design 
models are one step better in their simulation of the entire engine and the designer can use them with 
more confidence. After discussing with engineers at Stirling Technologies Co. (now Infinia, our teaming 
partner and one of several companies with which we communicate), we chose to focus on the flow exiting 
the heater channels and entering the expansion space, as it turns radially inward, then reverses. 

This work provides a complete description of oscillatory flow in a geometry that is representative of a 
Stirling engine expansion space. Secondarily, these investigations generated a data set that could be used 
as a standard for Stirling engine CFD simulations. It is a fact that oscillatory fluid dynamics data of this 
nature are generally lacking and hence there is a need for this type of investigation. An exception is an 
MIT study by Kornhauser 4 that presents wall heat flux measurements in a heat exchanger/two-space 
experiment. Further exceptions include work carried out at the University of Minnesota in the past by 
Seume and Simon, ' ’ Seume et ah, Friedman and Simon," Seume et al. and Qiu and Simon. ' The 
most recent work was, and is currently being carried out by Niu et ah 5 ' 6 and Simon et ah 14 This work 
focuses on heat transfer in a Stirling engine regenerator. 

To this end, this section presents hot-wire anemometry and flow visualization data from an 
experiment that replicates important features of oscillatory flows in Stirling engines. These features 
include impinging and sink flows, large-scale separation zones, recirculation bubbles, unsteady shear 
layers and temporally varying transition and relaminarization zones. The test operating conditions (see 
below) were chosen to match representative Reynolds and Valensi numbers by scaling the test section up 
and the oscillation frequency down. 


2.2 Nomenclature 


Symbol 

Description 

Units 

BDC 

Bottom dead center 

[no units] 

CAD 

Crank angle degrees 

[degrees] 

f 

Data sampling frequency 

[1/s] 

r 

Radial variable 

[m] 

R e max 

Maximum Reynolds number 

[dimensionless] 

RTD 

Resistance temperature detector 

[no units] 

RPM 

Crankshaft revolutions per minute 

[1 /min] 

TDC 

Top dead center 

[no units] 

D 

Disc diameter 

[m] 

S 

Disc spacing 

[m] 

TI 

Turbulence intensity 

[dimensionless] 

II u II 

Velocity magnitude 

[m/s] 

u 

Vector velocity 

[m/s] 

Va 

Valenci number 

[dimensionless] 

X 

Axial variable 

[m] 

e 

Crank angle 

[degrees] 

CO 

Crankshaft rotational speed 

[1 /min] 
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2.3 Experimental Set-Up 


The test facility is shown schematically in fig. 2.1. A Scotch yoke drive is used to develop sinusoidal 
piston motion of a constant rotational speed, co, where oo is variable from 4 to 120 revolutions per minute 
(rpm). Attached to the drive are a top dead center (TDC) indicator, a shaft angle encoder and a 
tachometer. The test section shown in fig. 2.1 is enlarged in fig. 2.2. The delivery tube inner diameter, Z), 
matches the test rig cylinder bore. The two large discs have a major diameter of 1 .22 m. The disc diameter 
is large to minimize the flow of momentum in and out of the test section. The delivery tube attached to 
the right disc in fig. 2.2 is 0.25 m in length. When the drive is at TDC, (the closest that the piston comes 
to the test section) the piston face is flush with the tube inlet plane (as shown in fig. 2.1). The piston 
stroke is 0.36 m (14 in.). The distance between the two large discs, S, is fixed for each case, but variable 
from 0.025 m (1 in.) to 0.127 m (5 in.). 


2.4 Data Collection 


2.4.1 Flow Visualization 

Pictures of the flow are captured with a 35 mm camera. The flow is seeded with particles, which, in 
turn, are illuminated with a laser light sheet. Particles are introduced into the flow with a Sage Action 
Model 33 Multi-Head Bubble Generator. This device combines compressed air, helium and a soap 
solution to produce nearly neutrally buoyant bubbles. The bubble size is adjustable from 6 mm to less 
than 1 mm. The bubbles used in the experiments are generally on the order of 3 mm. The bubbles follow 
the flow and give a 3-D, real-time view of the fluid motion. To illuminate a 2-D slice of the test section, a 
Melles Griot, He-Ne laser was used in conjunction with a cylindrical glass rod. When the beam passes 
through the rod, it spreads into a sheet of light, which, when appropriately placed, can illuminate any 
plane of interest within the test section. Attached to the end of the laser is an electronically controlled 
shutter which receives signals from the TDC indicator and the shaft angle encoder. The bubbles within 
the sheet are photographed with a Nikon FE2 35 mm camera loaded with ASA 400 or ASA 125 black and 
white film. 

The experimental procedure commences by selecting a disc spacing and an rpm. A crank angle during 
which flow is visualized is chosen, e.g., 15° to 30° (where 0° corresponds to TDC, 0° < 0 < 180° 
corresponds to the intake stroke, 180° corresponds to bottom dead center (BDC) and 180° < 0 < 360° 
corresponds to the blowing or exhaust stroke). The shutter controller is programmed to open the shutter at 
15° and close it at 30° (in the 15° to 30° example). This gives a 15° window in which the flowing 
particles within the sheet of light are illuminated. The shutter controller is programmed to open (and 
close) the shutter 100 times over this interval. The room is then blackened and the camera is set up to 
view the light sheet continuously (most manual cameras have a setting called ‘bulb’ which holds the 
camera shutter open for as long as desired). After 100 cycles, the camera shutter is closed and the result is 
an ensemble-averaged photograph of a 15° segment of crank angle duration. The process is repeated for 
other intervals in the cycle until the entire cycle is characterized. The flow visualization field of view is 
shown in fig. 2.3. Three different rpm and disc spacing operating points are documented, as summarized 
in table 2.1. 


TABLE 2.1.— CASES TESTED 


Case 

rpm 

Re 

lvt m;iv 

S 

(mm) 

Va 

Location 

SID 

I.a 

30 

7600 

54 

2300 

Tube 

0.251 

I.b 

30 

7600 

54 

2300 

Channel 

0.251 

Il.a 

30 

7600 

127 

2300 

Tube 

0.591 

Il.b 

30 

7600 

127 

2300 

Channel 

0.591 

III. a 

70 

17700 

54 

5400 

Tube 

0.251 

IILb 

70 

17700 

54 

5400 

Channel 

0.251 
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2.4.2 Hot-Wire Anemometry 

A TSI model 1053 hot-wire anemometer bridge was used for velocity measurements. The signal from 
the anemometer was sent to an IOTech A/D converter and then directed to a PC for recording. The 
angular position of the drive (from which piston position is calculated) was determined with an optically- 
triggered TDC indicator and a tachometer. Barometric pressure was recorded with an in-house pressure 
sensor. Ambient temperature was monitored with an RTD. The position of the anemometer probe was 
controlled with a Velmex stepper motor and an integrated slide unit. All data collection, probe positioning 
and data reduction were carried out with C code. 

Velocity and turbulence intensity (77) profiles were documented by traversing the probe within the 
test section. A single traverse was taken within the delivery tube and several traverses were taken within 
the disc space (see fig. 2.4) 

Disc space traverses are (following the convention of fig. 2.4) in the x-direction. The single traverse 
in the tube space was in the /'-direction. The coordinate origin is taken on the tube centerline (the “Axis of 
Symmetry”) at the target (left in fig. 2.4) disc. For reference, the disc on the right in figure 2.4 will be 
named the ‘back’ disc. The sequence of spatial points in a given traverse was randomized to prevent any 
temporal changes from being interpreted as spatial trends. 

At any given point, the analog signals from the anemometer and the TDC indicator were sampled at a 
frequency,/ of 500 Hz. Given/, ot> and the signal from the TDC indicator, it was possible to determine 
the location of the crank and piston. The raw data were ensemble averaged over 100 cycles to produce 
velocity and 77 values resolved in time within the cycle at each measurement station. For example, in the 
30 rpm case, the total run time for a single point was 200 seconds (100 cycles). Thus the total number of 
samples was 10,000. These 10,000 points were then separated by crank angle (the resolution was 1/2 
CAD) and the result was averaged over the 100 cycles. The probe was then moved to a new station and 
the process was repeated until an entire profile was completed. The cases investigated are shown in 
table 2.1. 


2.5 Results 

Unlike section 1 of this work, under oscillatory flow, it isn’t so simple to compartmentalize results 
into different subheadings (i.e., flow visualization, average velocity profiles, etc.) and discuss each 
separately. For example, flow visualization results give qualitative insight into the fluid behavior within 
the test section. Of course, little is revealed in terms of velocity magnitude aside from whether the flow is 
‘fast’ or ‘slow.’ This information is retrieved with hot wire velocity data. A compilation of hotwire data 
produces velocity and 77 plots. However, as we will shortly see, there are regions where 77 values not 
only exceed 30% but sometimes are as high as 100% or more (The magnitude of 77 values and their 
implications will be discussed shortly). We then turn to visualization and say; “of course, the wire is in a 
separation zone, and in these regions there is no reason to believe that an average velocity as indicated by 
a hot wire reflects anything about a typical flow condition.” Thus, in this section we present and discuss 
the results as a collection of flow visualization data, ensemble averaged velocity data and turbulence 
intensity data. 

The three cases of table 2.1, each with a sub-case which presents data in the tube and a sub-case in the 
disc space, were conducted. All six of the entries in table 2.1 will be described below. 

To the reader: you may find the following discussion conspicuously scant of flow visualization 
results. This is intentional; 144 pictures were taken to cover the three cases. A complete presentation 
would necessarily be lengthy and verbose. Thus, we present summary figures (sketches of the major 
features) of the results. However, we wouldn’t want to disappoint; to this end, we present a subset of the 
flow visualization photos in section 4 - Comparisons Between Experimentally Gathered Results and CFD 
Generated Data. We feel that, when viewed in light of CFD results, the flow visualization results ‘speak 
for themselves.’ 
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Please note that a hot wire is unable to resolve flow direction, hence the quantity 


llu II 



appears in velocity plots (where u, v and w are the vector velocity components). Suppose that the wire is 
in a flow and the flow is moving with some vector velocity u . Suppose some time later the flow reverses 
and has a vector velocity U that is anti-parallel to u with the same vector magnitude. The anemometer 
will report identical velocities. The point is that to a hot wire, which essentially measures the rate of 
cooling of a small cylinder, positive and negative velocities look the same. Thus, to get a grasp of the 
flow direction, one must look at the corresponding crank angle. 

We wish to acquaint the reader with the type of hot wire data to be presented. First of all, we are 
investigating average velocities, II u II, or turbulence intensities, 77, 



as a function of both crank angle, 0, normalized axial location x/S and normalized radial location, r!D. In 
this light, we could present a collection of two-dimensional slices of the data set with comments about the 
connection between each of the slices. Owing to the fact that crank angle resolution is 1/2 of a degree, this 
method would lead to an excessive number of figures. 

In this light, we choose to present three-dimensional plots with the independent variables presented as 
functions of two of the three independent variables, with the third independent variable held fixed. We 
begin by referring the reader to figure 2.5: this is the 3-D view that Matlab produces when it is fed a 3-D 
data set; ||u|| vs. x/S with r/D held constant in this case. Figures 2.6 and 2.7 are pictures of the same data 
set from different angles. The point is that, although a complete picture is presented in one figure, the 
reader has to imagine how it may look from different angles. 

Another note; the variable 0 appears in the discussion below. The reader should be aware that 0° < 0 
< 180° corresponds to the “suction” stroke, i.e., the piston is drawing fluid into the test section. On the 
other hand, 180° < 0 < 360° corresponds to the blowing half-cycle. In the following, we use the words 
blowing, blowing half-cycle, suction, suction stroke, etc. to represent the above stated crank angle 
regimes. Finally, we give the crank angle 0 = 0° = 360° the special name ‘Top Dead Center’ or TDC, and 
the angle 0 = 180° the name ‘Bottom Dead Center’ or BDC. 

One further note about average velocities II u II and turbulence intensities 77; without exception, in the 
77 profiles that will be presented shortly, there are regions where TI exceeds 30%. According to the 
literature, (Smits, 15 Hinze 3 ), when 77 values exceed 30% there is instantaneous flow reversal across the 
hot wire. Recall that the wire acts as a rectifier, thus, the net effect is that in zones of high values of 77 
measured average velocity values are too high and measured TI values are too low. For example, 
suppose that two velocities were taken and the actual value of velocity was 3 . 1 m/s and -3 m/s. The 
actual average velocity over these two points is 0.05 m/s and the actual TI over these two points is 
8630% (this ridiculous TI value is strictly for demonstration purposes). The anemometer will report 
velocity values of 3.1 and 3 m/s. Therefore, the average velocity II u II, according to the anemometer, is 
3.05 m/s and according to the anemometer, the TI is 2.3%. 

According to Smits, 15 for measured TI values of 30%, the error associated with TI measurements is 
10%, and the error associated with average velocity measurements is 10%. For measured TI values 
greater than 30%, the average velocity measurement error grows strongly with TI whereas for 77 less than 
30%, the errors in TI and average velocity measurements are small. 
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Case la— Tube flow, Va = 2300, S/D = 0.245 

We begin the explanation with a presentation of ensemble average velocity and 77 data. Shown in 
figure 2.8 is average velocity plotted against crank angle and normalized radial location. One sees that the 
flow in the tube follows the sinusoidal motion of the piston. The changing boundary layer thickness can 
be observed near r/D = -0.5 and some streaking of the flow is seen near the centerline ( r/D ~ 0) during the 
exhaust stroke (180 to 360 CAD). Presented in figure 2.9 is a 77 plot. With the exception of BDC (and 
the surrounding 0 values) 77 values are low (77 < 0.3) and we may assume that velocity measurements are 
accurate. 

Turning to flow visualization, we present a summary report of crank angles near TDC on the intake 
stroke in figure 2.10. In the early portion of the cycle, the flow has little momentum, and is reversing. As 
the cycle proceeds, flow is drawn into the cylinder in an apparently laminar manner. One curious and 
noteworthy flow feature during this portion of the cycle is that visualization reveals no growth or presence 
of separation bubbles (fig. 2.1 1) for this portion of the cycle. At 45 CAD, separation bubbles appear in the 
visualization results and at about 80 CAD, the flow seems to trip to turbulence (see fig. 2.12). An 
interesting feature of the flow is that for 45° < 0 < 105° there is a separation bubble that does not grow 
with increasing crank angle (i.e., stronger suction or pressure gradient). For illustrations of these features 
(see figs. 2.12 and 2.13). As the suction stroke proceeds, the flow remains turbulent and the separation 
bubbles begin to grow (fig. 2.14). At BDC, the piston achieves zero velocity and consequently, the flow 
in the tube becomes full of small-scale turbulence. Upon reversal, the flow begins as a turbulent slug, and 
is blown out of the cylinder (fig. 2.15). Following in the wake of the turbulence is a laminar-like flow. 

The turbulence has apparently been dissipated by the walls and the strong temporal acceleration. As the 
piston continues the exhaust stroke, all of the turbulent flow is driven out and a laminar-like regime is 
maintained almost until TDC, as shown in figure 2.16. 

Case lb — Disc space flow, Va = 2300, S/D = 0.245 

This is a collection of measured velocity and 77 profiles and flow visualization results in the disc 
space. We begin with a presentation of flow visualization results. As a first comment, we note that during 
the blowing stroke, close to both TDC and BDC, the flow is more or less quiet in nature (i.e., lacks any 
perceptible momentum). Flowever, for intermediate crank angles, the flow is drawn radially inward in a 
more or less uniform fashion. This idea is presented pictorially in figure 2.17. The blowing stroke is 
different in nature; the flow feels the effect of the piston-driven pressure gradient almost immediately. A 
jet quickly forms on the target wall and with it comes a separation/recirculation bubble at the sharp edge 
of the tube exit (see fig. 2.18). The jet undergoes transition and radially outward of the transition point is 
a turbulent zone. The jet loses it’s coherent nature at a normalized radial location of approximately r/D = 
1.2. These comments are illustrated in figure 2.18. 

Turning to anemometry results, velocity profiles are presented at normalized radial locations, r/D, of 
0.68, 0.81, 0.93, 1.05, 1.17, and 1.30 in figures 2.19 through 2.24. One can see that for smaller values of 
r/D, during the blowing portion of the cycle, there is a coherent jet along the target wall (i.e., for x/S < 

0.4). As one moves radially outward, the jet thickness first decreases, then begins to increase. For r/D 
greater than 1.17, it appears that the jet spreads out to fill the entire channel (see figs. 2.23 and 2.24). 

Looking at 77 profiles (figs. 2.25 through 2.30) we see that, indeed, there is a region of low turbulence 
in the wall jet region until r/D = 0.93. For r/D greater than 1.05, the turbulence intensity during the 
blowing cycle is approximately uniform across the center of the channel. However, the 77 near the walls 
is much larger. This suggests a thick boundary layer as a result of the strong temporal deceleration of the 
flow. The 77 levels indicate a region of energized, turbulent flow. 

During the suction stroke (see figs. 2.19 through 2.30), the velocity profiles exhibit a uniform sink 
flow appearance. Moreover, for all radial locations, except r/D = 0.93, we have the fortunate condition 
that away from TDC (perhaps at 75 CAD or so) 77 levels are of a reasonable value. It is likely that in the 
late portion of the blowing stroke, a fair amount of residual turbulence is deposited in the disc space, 
which remains during the initial phase of the intake stroke. As the stroke proceeds, perhaps the turbulent 
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residuals from the previous half-cycle are re-ingested in the tube, and, as the flow accelerates, it stabilizes 
(i.e., becomes laminar-like). In any case, the combination of flow visualization and low TI values seems 
to suggest that velocity profiles in this regime are trustworthy. 

Case Ha— Tube flow, Va = 2300, S/D = 0.577 

Figure 2.31 shows that during the early part of the intake portion of the cycle, 0 < 0 < 180, the 
velocity near the tube walls is of a higher absolute value than the core flow velocity. This suggests that 
the higher momentum fluid in the core flow (see values for 0 < 270°) is slower to respond to change in 
pressure gradient imposed by the piston. Viscous forces near the wall reduce the momentum of the near- 
wall fluid and allow it to respond more quickly to the change in pressure field. 

The exhaust portion of the cycle 180 < 0 < 360 is somewhat different. Here, the flow is being driven 
out of the tube by the piston. It seems that for this case, the flow behaves like a developing pipe flow, 
with the core flow velocity being flat and following the piston velocity. Boundary layer growth can be 
seen, particularly for the deceleration portion of the cycle. 

Figure 2.32 presents a plot of turbulence intensity in the tube space. Periods of strong acceleration 
(e.g., 0° < 0 < 45°) show dissipation of turbulence, whereas periods of deceleration (e.g., 90° < 0 < 135°) 
show large production of turbulence. One may also notice that, excluding 160° < 0 < 210°, the absolute 
value of TI is low; hence we may consider average velocity data outside of this range to be reliable. 

We now turn to flow visualization. Figure 2.33 shows that near TDC on the intake stroke, the flow 
has a low over-all momentum and is locally and instantaneously reversing. As the cycle proceeds, flow is 
drawn into the cylinder and, in the early portion of the cycle, seems to exhibit a ‘alternating’ behavior, 
perhaps from cycle to cycle, with respect to the tube centerline. This is to say that on some strokes the 
flow seems to be drawn from the top of the disc space whereas on other strokes it seems to come from the 
bottom. Also noteworthy is the presence and growth of separation bubbles at the sharp-edged tube 
entrance. These phenomena are sketched in figure 2.34. As the cycle proceeds, the alternating behavior 
seems to disappear. In the degree range 60° < 0 < 90°, the flow seems to be axisymmetric, with a laminar 
core and growing separation bubbles (fig. 2.35). After 90 CAD (maximum intake piston speed) the flow 
appears to rapidly trip to turbulence. Again, separation bubble growth is present, but the core flow is 
turbulent in nature (fig. 2.36). As the intake stroke proceeds to completion, the flow decelerates with the 
piston, and reverses. The early part of the blowing stroke shows that the bulk flow is fully turbulent 
(fig. 2.37) and remains so until 240 CAD. After 240 CAD, the flow seems to Taminarize’ and it remains 
stable as shown in figure 2.38. This is perhaps due to the sustained exposure to temporal acceleration and 
dampening by the tube and piston walls. At the end of the blowing stroke (not shown), flow deceleration 
causes turbulent structures to form in the tube flow which persist into the next cycle. 

Case lib — Disc space flow, Va = 2300, S/D = 0.577 

The area of concern is the disc space (see figs. 2.3 and 2.4). Flow visualization reveals that during the 
blowing stroke and a portion of the suction stroke, there is a large standing vortex (or eddy) between the 
discs. The blowing stroke energizes the eddy; the suction stroke simply moves the eddy slightly. The 
motion of the eddy, along with other pertinent flow features is sketched in figures 2.39 and 2.40. 

Figure 2.39 shows that during the blowing stroke there is a strong and coherent wall jet. This jet generates 
an eddy, which initially forms at the sharp edge of the tube exit and, as crank angle increases, its center 
follows the dashed red line upwards. The region in the wake of the eddy and to the left of it is more or 
less dead (motionless) flow. Figure 2.40 is a summary of flow visualization results for the suction stroke. 
The eddy from the blowing stroke is pulled downwards in the general direction of the bulk flow but 
rapidly loses its momentum and becomes very weak after about 30 degrees of crank angle. The main flow 
which fills the volume that is opened is pulled from the back wall where the momentum from the 
previous half-cycle is low, with velocity seeming to decrease as one moves rightwards. 

Turning to anemometry data, we present a collection of profiles in the disc space at normalized radial 
locations, r/D, of 0.68, 0.81, 0.93, 1.05, 1.17, and 1.30. As a first comment, for all profiles with r/D > 
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0.81, TI never drops below 30% (with the exception of some regions at the radial location r/D = 1 .30). In 
short, we must look on ensemble-average velocity profiles with a somewhat skeptical eye. Turning to 
figure 2.41 we present a velocity profile at r/D of 1.30. At this radial location, there is some remnant of 
the wall jet at 0 = 0° and x/S ~ 0.10, but more or less, velocities are low, turbulence intensity is high (see 
fig. 2.42), and for all practical purposes measurements of velocity in this flow are too high in absolute 
value. We move radially inward, and present velocity and TYprofiles at r/D of 1.17 (figs. 2.43 and 2.44). 
Again, we see a similar behavior; remnants of a wall jet are present, and again TI is high and velocities 
are low. Next we present profiles at r/D of 1.05 (figs. 2.45 and 2.46). Here we have the fortunate 
condition of having flow visualization results to aid in the discussion. As is the case for larger radial 
locations, we still have the condition that TI is too high for accurate anemometry; the recorded velocities 
are too high and the TI values are too low. However, the anemometry data show the presence of a wall jet. 
The jet thickness is approximately 30 to 40% of the channel width. Looking at figure 2.47, we see that 
there is, in fact, a wall jet. We can expect that it is as thick as the hot wire data says. Turning back to our 
summary figures (figs. 2.39 and 2.40) we may expect a wall jet throughout the blowing stroke and a 
relatively weak sink flow on the intake stroke in the vicinity of the back wall (i.e., for large values of x/S). 
In figure 2.45, the jet is visible. It even seems to maintain some of its momentum throughout the first 30 
degrees of the suction stroke. Turning now to the idea of a sink flow in the vicinity of x/S =1, there exists 
a region of high velocity fluid in this area, which suggests that this flow may be present. We continue 
inward now and present profiles for r/D of 0.93 in figures 2.48 and 2.49. Again, the same behavior is 
observed as for the larger radius cases. Moving inward again, we present profiles at r/D = 0.81 in figures 
2.50 and 2.5 1 . For the first time, we see that TI levels drop below 20%. This happens near x/S = 1 and x/S 
= 0, during periods of high flow momentum (see the top left and bottom right portions of fig. 2.5 1 ). For 
these data the hot wire values are accurate. Moving in again, we present profiles for r/D = 0.68 in figures 
2.52 and 2.53. We also present flow visualization results in figure 2.54. We see a coherent (thicker, 
approximately 60% of the channel width) wall jet and a definite sink flow near x/S = 1 . 

Case Ilia— Tube flow, Va = 5400, S/D = 0.245 

The field of view is the tube portion of the test section (see fig. 2.3). At the beginning of the cycle ~0° 
—>18° after TDC the piston is drawing fluid into the field of view from the right. In this degree range, it 
is observed that the fluid acceleration is spatially non-uniform (see fig. 2.55). After this degree range, 
from ~18° — > 72° after TDC, a low-turbulence flow forms in the core of the pipe and separation bubbles 
grow downstream of the sharp inlet. Figure 2.56 shows a sketch of the separation zone on the right and 
how it grows in this degree range. Continuing through the cycle, a different behavior is observed from 
~72° — > 126° after TDC. The flow turbulence remains low and the separation bubbles continue to grow, 
but an alternating flow behavior is observed. To clarify, it appears in the photos that on some cycles, flow 
is drawn from the upper part of the disc space whereas on other cycles it is drawn from the lower part. 
This leads to a correspondingly variable separation zone from cycle to cycle. A sketch of the general 
features observed in this degree range is shown in figure 2.57. At this point, it is noted that for all values 
of crank angle between 90°-> BDC, the flow is decelerating. The deceleration destabilizes the flow and 
produces large-scale turbulence. At -126°, large scale and strongly turbulent activity becomes prominent. 
With this, the growth of the separation bubbles changes. The growth rate in the streamwise direction of 
the bubbles on both top and bottom surfaces remains relatively constant. However, the bubble on the 
upper surface grows rapidly in the cross-stream direction (see fig. 2.58). The bubbles grow until BDC is 
reached. At BDC, the tube is entirely filled with eddys of multiple scales and the bulk convection is veiy 
low. After BDC, the drive pushes the piston and expels fluid from the cylinder to the disc space. The 
pictures of the exhaust stroke, e.g., BDC — > TDC, show a “wave” propagating from left to right 
(fig. 2.59). This process may be understood if one assumes that there are two “slugs” of fluid within the 
test section during the exhaust stroke. First, there is a low-turbulence slug of fluid that is close to the 
piston face. This low-turbulence slug is trapped within the tube and simply shuttles back and forth. 
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Second, there is a slug of fluid that escapes into the area between the two discs and is shuttled between 
this space and the cylinder. The above idea is illustrated in figure 2.59. The second slug experiences flow 
separation and is highly turbulent in nature. Apparently, the turbulence in this slug is not effective in 
agitating the near-piston slug. No separation bubbles form within the cylinder during this part of the 
cycle. A perplexing flow feature is observed during the exhaust stroke at the very end of the stroke. From 
-162° ABCD — > TDC° there is a non-uniform velocity profile within the tube. In the upper left-hand side 
of the tube, the flow has a higher velocity. This corresponds to the entry flow at the beginning of the next 
cycle. This is sketched in figure 2.60. 

A velocity profile is presented in figure 2.61. We see the flow follows the piston in a sinusoidal 
fashion. One immediately notices that the profiles on the intake and exhaust strokes have different shapes; 
the profile on the exhaust stroke is ‘flat’ whereas the profile on the intake stroke is convex. One may 
attribute the difference in shape to the fact that on the intake stroke, separation bubbles form at the tube 
entrance; hence the flow must accelerate around them. On the exhaust stroke there is no separation bubble 
influence, hence the flow assumes a flat profile and follows the piston. Figure 2.62 is a plot of turbulence 
intensity; as one may expect, periods of strong acceleration (25° - 90° and 200° - 270°) laminarize the 
flow. The careful reader may notice that on the intake stroke, turbulence quantities begin to increase 
around 120 CAD, whereas, on the exhaust stroke, a similar behavior is not seen until 350 CAD. In other 
words, 77 values do not follow a sinusoid. This may be attributed to the fact that on the intake stroke, 
fluid is being drawn in from the disc space, which is an inherently turbulent flow. 

Case Illb — Disc space flow, Va = 5400, S/D = 0.245 

We begin with flow visualization results. In the early portion of the suction stroke (0° < 9 < 18°) the 
flow in the disc space seems to lack any coherent pattern; it is more or less a field of small-scale turbulent 
eddies (see fig. 2.63). As the piston continues the suction stroke, fluid begins to move radially inward. For 
small crank angles, fluid near the sharp-edged tube entrance begins to flow towards the tube (as shown in 
fig. 2.63). One can see that as crank angle increases, the pressure gradient is felt quickly in the cross- 
stream direction (blue arrows for small crank angles and red and blue arrows for larger ones). A picture of 
larger crank angles is shown in figure 2.64. Somewhere between 60 and 70 CAD, the sink flow reaches 
the target wall. In this region, as far as can be told, the velocity profile is more or less uniform. At the 
same time, the effects of the pressure gradient are being felt at larger radial locations. By the time 
maximum piston speed is achieved (i.e., 90 CAD), the flow is approximately uniform across the channel 
and covering the entire span of the flow visualization field of view (fig. 2.65). For the duration of the 
suction stroke, the flow remains uniform throughout the channel with flow velocity and piston velocities 
decreasing. 

At BDC, the flow is, as is reasonable, devoid of momentum. As the piston begins to blow fluid out of 
the cylinder, a wall jet and separation bubble rapidly form. These features, and their development, are 
shown in figure 2.66. Everywhere adjacent to the separation bubble and near the target wall, the fluid has 
a remarkably high momentum as it streaks radially outward. The bubble, at its thinnest cross-stream 
dimension, occupies one-half of the disc space; at its thickest point, it occupies approximately seven- 
eighths of the channel. Figure 2.66 shows the veiy end of the separation bubble growth, (i.e., 0 ~ 310°). 

At smaller crank angles of the 190° < 0 < 310° range, the bubble is smaller, and the ‘streaks’ in the 
channel near the target wall are at a smaller radial location. For crank angles greater than 310°, the flow 
destabilizes, and except very near the tube exit, the channel fills with turbulence (fig. 2.67). 

We now turn to anemometry data, and as the reader may have remarked while looking at the flow 
visualization summary diagrams, all of the hot wire sampling stations are within the flow visualization 
field of view. We begin the presentation of anemometry results by discussing a velocity profile at r/D = 
0.68 (fig. 2.68). We see that during the suction stroke, as the flow visualization indicated, the flow 
responds quickly near the tube exit (near x/S = 1). Shortly thereafter, the flow takes on a uniform sink 
pattern. The flow follows the piston motion and slows in the vicinity of 180° and, then, during the 
blowing stroke, a strong wall jet forms and lasts for the rest of the cycle. Turning to 7/ data (fig. 2.69), we 
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see that there are two regions where TI values are at or below 30%. These regions correspond to the wall 
jet, shortly after the piston has reached its maximum exhaust velocity and the channel-spanning sink flow 
shortly after the piston has achieved its maximum intake velocity. This corresponds nicely with figures 
2.65 and 2.66. 

Moving radially outward, we present a velocity plot at rID of 0.8 1 in figure 2.70. We observe a 
similar behavior to the smaller radial location; there is an excess of momentum in the degree range 270° < 
0 < 360° and 0° < 0 < 100° near the back wall (x/S ~ 1). One might suspect that the separation bubble that 
forms on the back wall is influencing the hot wire signal, and, though the signal indicates a relatively high 
velocity, it may be due to large turbulent eddy transport. To clarify this idea we turn to figure 2.71. We 
see that the TI values in these regions are relatively high (30 to 40%); however, this does not clarify 
weather there is a high-energy eddy present there or if the flow is in fact streaking along the back wall. 

Another step out lands us at r/D = 0.93. A velocity profile is presented in figure 2.72. Again, a wall 
jet and uniform sink flow are present. Also present is the region of high indicated velocity near the back 
wall. Turning to figure 2.73, we find uniform and high levels of TI across the channel at 0 = 300°. The 
flow visualization does not indicate the presence of a strong flow in these regions. It indicates the 
presence of a separation bubble (or perhaps its remnants). Hence, we attribute these high-indicated 
velocities to a recirculating flow at this time; outward near the target wall, and inward near the back wall. 

The same trends that were discussed above persist as one continues radially outwards. The results are 
presented in figures 2.74 through 2.79. One thing to notice in these figures is that at larger radial 
locations, the flow is more turbulent in nature. In general, the regions that are suitable for hot wire 
anemometry remain stationary (i.e., in the same crank angle and spatial location) but decrease in size with 
larger radii. One anomaly is present in the larger radial location plots; in the vicinity of 240° for 0.93 < 
r/D < 1 .03, a low-turbulence region forms near the target wall and, with an increase in radius, grows 
towards the back wall. This is during the blowing stroke as the piston is approaching its maximum 
velocity. It corresponds to a ‘ridge’ forming on the velocity profiles (e.g., see fig. 2.74). Perhaps this 
‘ridge’ corresponds to the flow accelerating around the separation bubble. 

2.6 Wall Shear Stress as a Function of Crank Angle 

Of interest was the determination of wall shear stress as a function of crank angle. This information is 
useful in the sense that it is the next level of complexity in code validation. In this section, we present 
wall shear stress measurements, as derived from the ensemble average velocity profiles presented above, 
as a function of crank angle. We will also outline the difficulties in measuring shear stress with a hot- 
wire, and the errors associated with the measurement. 

To begin, we refer the reader to figures 2.99 and 2.100. These figures show that there is an effect on 
the hot-wire as the probe approaches the wall. This is unfortunate in the sense that wall shear is described 
as: 


X wall ~ 


du 

dx 


wall 


Thus, to directly evaluate wall shear stress, one needs to know the gradient at the wall, which is 
influenced by the wall effect. To circumvent this somewhat inconvenient problem, we turn to theory. 
According to Burmeister, 1 in a turbulent flow, near a solid boundary, the following equation (a.k.a. the 
law of the wall) holds: 


+ + 
a = x 


where 
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Thus, when wall shear stress is known, and the velocity profde near the wall is cast in these non- 
dimensional variables, the velocity profde near the wall is linear with an intercept at the origin. Since we 
want to determine wall shear stress, this presents a small problem. To this end, we guess a value of wall 
shear stress and then check this guessed value against the collected data. We now present the algorithm 
used to accomplish this feat and refer the reader to appendix C for the Matlab code used to cany out the 
computations. 

2.6.1 An Algorithm Capable of Determining Wall Shear Stress from Corrupted Experimental Data 

Step 1. Truncate each velocity profde (at each crank angle) to eliminate far- wall velocity data. The 
truncation is made at x/S = 0.1 for all cases and crank angles. We feel that at larger axial locations (larger 
x/S) the “law of the wall” as described by Burmeister, 1 does not hold. 

Step 2. With the remaining data, 0 < x/S <0.1, locate where the wall influence on the hot-wire 
disappears as one traverses away from the wall. With reference to figure 2.100 one can see that the 
profiles in the near-wall region (x/S < 0.008 or so) are corrupted regardless of crank angle. The simplest 
way to remove the wall effect is to look for the minimum value of velocity at each crank angle in the 
profde and remove all points inbord of this velocity minimum. 

Step 3. Having now eliminated far-wall data and near wall effects, we are ready to guess a slope and 
hence a wall shear stress. This is achieved by first assigning a fixed point at the origin. We let velocity 
equal zero at x = x/S = 0. Next we take the minimum value of velocity (as identified in step 2) and the 
next five points outboard from the minimum. We then take their average in both velocity and position. 
The result is a single point in space and velocity corresponding to the average value of the near- wall 
profde. 

Step 4. Fit a line through the origin and the average value determined in step 3. Multiply the slope of 
this line by viscosity to get shear stress. 

We now have to identify any errors associated with the truncation and averaging of the profiles to 
determine shear stress. 

Step 5. Now that shear stress is known from Step 4, we take the minimum value of velocity, and the 5 
nearest points outboard of it (from Step 3), and convert them to u + , x + coordinates. 

Step 6. For each matched u + , x + we define the quantity: 


RSS = 



and use this as an error bound on the slope (and hence wall shear stress) derived in Step 4. 
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2.6.2 Wall Shear Stress for Case lib (Re max = 7600, Va = 2300, S = 127 mm, S/D = 0.591) 

Given average velocity data for Case lib and the algorithm above, we present the results in 
figures 2.80 through 2.85. For r/D = 0.68 (fig. 2.80), we see that the wall shear stress is low for a large 
portion of the intake stroke (0 < 0 < 180). The large stresses in the early portion of the intake stroke are 
likely due to the fact that in this case, there is a strong remnant of the wall jet present from the previous 
half-cycle. During the blowing stroke, (180 < 0 < 360), there is a more or less steady rise in shear until 
about 310 CAD. At this point the shear increases dramatically and might indicate a transition to turbulent 
flow. Figure 2.52 shows a rapid change in the wall layer velocity profile at 0 = 300°. We see that for this 
particular radial location the error, as defined above, is generally in the neighborhood of 4-9%. Moving 
radially outward, in the next figure (fig. 2.81), we see a somewhat similar behavior at the radial location 
r/D = 0.81. Though difficult to know precisely, transition may be taking place when 0 = 300°. Figure 2.50 
shows a much thicker near wall layer at this radial position but a rapid change in shape at 0 = 300°. 
Another step out takes us to r/D = 0.93 (fig. 2.82). There are two things to note about this particular 
profile; first, there may be an indication of a short-lived transition at about 300 CAD. Second, the 
magnitude of the maximum shear stress has dropped by about a factor of three. Figure 2.48 shows that the 
velocities are lower at this radial position and, though there is a change in shape at 0 = 300°, the effect is 
weak. The next profile out (r/D = 1.05, fig. 2.83) is somewhat remarkable in terms of its ‘smoothness’ 
and that the shear stress values are quite low compared to the previous stations’ values. The flow has 
slowed considerably at this radial location. Figure 2.45 shows very low velocities and a thick near-wall 
layer. The next radial location (r/D = 1.17, fig. 2.84) again becomes jagged and transition is difficult to 
detect. Figure 2.43 shows a very low and flat near-wall velocity profile around 0 = 180°. The shear stress 
at this time may have been corrupted by occasional reverse flow. The last profile at r/D = 1.30 (fig. 2.85) 
is somewhat bizarre in its own right, and little can be said about it except that the value of shear stress is 
quite low. The oddity of the profile is not entirely unexpected, considering the nature of the average 
velocity profile at this point (fig. 2.41). Average velocities are quite low. 

2.6.3 Wall Shear Stress for Case lb (Re max = 7600, Va = 2300, V = 54 mm, S/D = 0.251) 

Using the algorithm described above, plots were made for the radial locations r/D = 0.68, 0.81, 0.93, 
1.05, 1.17, and 1.30 and are presented in figures 2.86 through 2.91. One thing to notice for the smallest 
radial locations (figs. 2.86 through 2.88), is that there is a distinct difference between the blowing and 
suction strokes. Namely, during the blowing stroke, 180 < 0 < 360, there is a considerably higher shear 
stress than on the suction stroke. This is in line with the hot wire and flow visualization results presented 
earlier. For these small radial locations (r/D = 0.68, 0.81, and 0.93), there is a strong wall jet during the 
blowing stroke, and on the intake stroke, flow is drawn from the entire channel, but more strongly from 
the area near the back wall. Thus, the shear stress is actually reversing during the cycle, though the hot 
wire measurements cannot resolve it. These characteristics would suggest shear stress plots like the ones 
presented in figures 2.86 through 2.88. One thing to note about these profiles (figs. 2.86 through 2.88) is 
that, with the exception of r/D = 0.93, the error is somewhat high on the blowing stroke. Hence, it is 
likely that general trends are being correctly described whereas numerical values of stress might be 
incorrect. For the larger radial locations (r/D = 1.05, 1.17, and 1 .30, figs. 2.89 through 2.91) the behavior 
is somewhat different. The magnitude of shear stress is approximately equal on the blowing and suction 
strokes. This is reasonable in light of the flow visualization and anemometry results discussed above; the 
best example is figure 2.24. One can see that the velocity profile is similar on both the intake and exhaust 
strokes - there isn’t a strong wall jet present on the blowing stroke. In these larger radius profiles, the 
error is quite low and we expect that the values of shear stress are correspondingly accurate. 

2.6.4 Wall Shear Stress for Case Illb (Re max = 17700, Va = 5400, V = 54 mm, S/D = 0.251) 

Presented in figures 2.92 through 2.97 are walls shear stress profiles at the radial locations r/D = 0.68, 
0.81, 0.93, 1 .05, 1.17, and 1.30. These profiles are somewhat similar to those presented above for Case la. 
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The value of shear stress is higher on the blowing stroke than it is on the suction for the closest four radial 
locations (figs. 2.92 through 2.95). However, these plots are somewhat fraught with error. More than 
likely it is the case that the increase in wall jet velocities has thinned the boundary layer and made 
accurate shear stress measurements difficult or impossible. For the last two profiles this is not the case 
(figs. 2.96 and 2.97, r/D =1.17 and r/D = 1 .30), the error drops, presumably the boundary layer thickens 
and the values of shear stress are more reasonable. For these cases, the velocity profiles in figures 2.68 
through 2.79 can be referred to for further clarification of the near-wall flow. In this case, as in the last, 
the near-wall flow and the shear stress are reversing during the cycle. 

2.7 Conclusions 

A data set that is useful for CFD code validation has been generated. The geometry under 
investigation is spatially simple, yet it generates flow complex enough to challenge current available 
routines. The code developer will find the above ensemble average velocity data, 77 quantities, and wall 
shear stresses indispensable in proving that the code “works.” Moreover, from the design point of view, it 
has been observed that tube flows are generally sinusoidal in nature and relatively independent of engine 
parameters. It has been observed that flows in the disc space (which approximately models a heater head) 
generally show the presence of a wall jet (on the blowing stroke) and a sink flow (during suction). These 
features depend strongly on engine geometry (disc spacing) and Valensi number. 

In light of sections 1 and 2, the reader may ask what is the connection? Section 3 - A Comparison 
Between Oscillatory Flow and Unidirectional Flow serves to answer this. The main subject of concern is 
to determine whether or not it is reasonable to expect similar fluid dynamic behavior from steady 
unidirectional flow and oscillatory flow in the same geometry. The impetus is the fact that there is a large 
body of work available for impinging jets, whereas, there are few data for oscillatory flows in Stirling 
engines. This knowledge is useful for both the engine designer and code developer alike. 

A comment is in order regarding CFD development; the data provided above provide insight into 
what type of modeling may be necessary in different engine regions. For example, in the tube flow it was 
observed that by and large the flow was axisymmetric and sinusoidal. This type of flow could be well 
described with one-dimensional design codes given the proper choice of parameters (i.e., Reynolds 
number, Valensi number, etc.). Likewise, in the wall jet, given proper normalization and averaging 
procedure, one could accurately describe bulk fluid behavior. (Some sort of jet thickness as a function of 
engine parameters should be developed, which, in the authors’ opinion, there is insufficient data available 
for this purpose.) Regarding multi-dimensional modeling, one can expect that large-scale flow features 
could be well described with RANS type modeling. However, development of proper turbulence models 
will be necessary. With the wrong turbulence model, one cannot reasonably expect accurate calculations 
of wall shear stress values, turbulence quantities, shear layer shape, etc. This will be discussed in further 
detail in section 4 - Comparisons Between Experimentally Gathered Results and CFD Generated Data. 

2.8 Error Analysis 

As a first comment, hot wires have a certain calibration error associated with them. This source of 
error is handled in detail in appendix B. Outside of specific calibration errors, there are errors associated 
with measurements of the above-described flow. We now present the errors that we noticed. 

One potential error associated with anemometry in an oscillating flow is as follows; the probe 
introduces a disturbance to the flow. We introduce figure 2.98 to clarify our meaning. With reference to 
the ‘Side View’ in figure 2.98, one must imagine that during some portions of the oscillation cycle, flow 
is moving from left to right (as drawn) whereas at other times it is moving from right to left. The point is 
that the wake the probe creates will eventually be blown back over the probe and wire. Thus, the 
introduction of the probe into the flow might lead to higher 77 measurements than actually present in the 
undisturbed flow. Moreover, shortly after reversal, as the wake is passing over the probe, velocity 
measurements may be perturbed in that the wake will ‘shield’ the probe from the bulk flow velocity. 
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Unfortunately, we have no direct way to quantify the effect of these errors, we can only state that they 
may exist. One could get a feel for the magnitude of the error by conducting a flow visualization 
experiment with the probe in the flow. This is difficult in the sense that hot wire probes are remarkably 
fragile and introducing particles into the flow would for visualization purposes invariably break the probe. 

A second source of error arises when taking near-wall measurements. As the active part of the probe 
(the wire) is moved towards a solid boundary, conduction effects seem to become important 
(see fig. 2.99). In general, within a critical distance, ‘A,’ the wire thermally communicates with the wall 
via conduction through the air. The net result is that as more heat is removed from the wire (due to the 
presence of the wall) the measurement of apparent velocity increases. Of course, the flow must obey a no- 
slip boundary condition at the wall; hence velocity should go to zero as the probe approaches the wall. A 
sample data set is presented in figure 2.100. Figure 2.100 is a particular (i.e., A in fig. 2.99 varies with 
rpm) example of what one observes that near the wall (x/S is near zero). Regardless of crank angle all of 
the curves ‘follow’ this trend. At about x/S = 0.008, the curves begin to deviate, and, for larger axial 
positions, the curves are not correlated in any special way. Thus, we assume that the correlated portion of 
the curves is characteristic of wall conduction effects. To ‘correct’ this problem, we simply discard the 
points in the correlated region. Thus, we take A (as mentioned above), in this case as 0.008 and discard all 
points that are closer to the wall than A. In general, A ~ 0.008 for the 30 rpm cases and A ~ 0.004 for the 
70 rpm case. As a note, this correction has been made to all of the velocity and 77 plots presented in this 
report. 

A third source of error concerning the hot wire corresponds to pulling the probe out of the flow. 
Again, a picture speaks better than words, so we reference the reader to figure 2.101. The type of error 
illustrated in figure 2.101 happens only near x/S = 1. At x/S = 1, the wire is essentially at the lip of the 
access hole. Looking at raw data (fig. 2.102), one sees that the velocity near this axial location is too high 
(almost like a jump discontinuity), hence, data points near to x/S = 1 are treated as erroneous and not 
presented in the plots above. 

Outside of anemometry measurement errors, there are errors associated with probe positioning. We 
state them now. 

In the disc space for S = 54 mm, the error an the axial direction, x, can be described as; 

x/S = x/S±8x/S 


where 


8x/S = 0.002. 

In the radial direction, r/Z), the associated positioning error can be described as: 


r / D = r / D ±8r / D 


where 


5r/D = 0.005. 

In the disc space for S = 127 mm, the error an the axial direction, x, can be described as: 

x/S = x/S±8x/S 


where 


5x/S = 0.001. 

In the radial direction, r/D, the associated positioning error can be described as: 
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r/D = r/D±br/D 


where 

8r/D = 0.005. 

Finally, in the tube space the error in the radial direction, r, can be described as: 

r ID — r / D ±8r / D 

where 


8r ID = 0.0006. 


2.9 Figures 

Scotch Yoke Mechanism 



Figure 2.1. — The University of Minnesota Scotch Yoke Drive Facility and Test Section. 
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Figure 2.2. — The test section. 
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Figure 2.4. — Sampling locations and coordinate conventions. 



Figure 2.5. — A sample data set - the best view of small values of x/S. 
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Figure 2.6. — The same data set as figure 2.5 but viewed from a different 
angle - the best view of the effect of crank angle. 



Figure 2.7. — The same data set as figure 2.5 but viewed from yet another 
angle - the best view of the region near x/S = 1 . 
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Figure 2.8. — Ensemble averaged velocity profile for Case la (tube space 
flow, Va = 2300, S/D = 0.251). 



Figure 2.9. — 77 contour map for Case la (tube space flow, Va = 2300, 
S/D = 0.251). 
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Figure 2.10. — Flow visualization summary for the degree range 0° < 0 <15° for Case la 
(tube space How, Va = 2300, S/D = 0.251). 
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Figure 2.1 1 . — Flow visualization summary for 15 to 45 CAD for Case la (tube space 
flow, Va = 2300, S/D = 0.251). 
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Figure 2.12. — Flow visualization summary for 45 to 75 CAD for Case la (tube space 
flow, Va = 2300, S/D = 0.251). 
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Figure 2.13. — Flow visualization summary for 75 to 105 CAD for Case la (tube space 
flow, Va = 2300, S/D = 0.25 1). 
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Figure 2.14. — Flow visualization summary for 105 to 180 CAD for Case la (tube space 
flow, Va = 2300, S/D = 0.251). 
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Figure 2.15. — Flow visualization summary for 210 to 240 CAD for Case la (tube space 
flow, Va = 2300, S/D = 0.251). 
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Figure 2.16. Flow visualization summary for 240 to 330 CAD for Case la (tube space 
flow, Va = 2300, S/D = 0.25 1). 
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Figure 2.17. — A summary of flow visualization results for the degree range 30° < 9 < 
165° for Case lb (disc space flow, Va = 2300, S/D = 0.251). 
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Figure 2.18. — A summary of flow visualization results for the blowing stroke for Case 
lb (disc space flow, Va = 2300, S/D = 0.251). 



Figure 2.19. — Ensemble averaged velocity profile at r/D = 0.68 for Case 
lb (disc space flow, Va = 2300, S/D = 0.25 1 ). 
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Figure 2.20. — Ensemble averaged velocity profile at r/D = 0.81 for Case 
lb (disc space flow, Va = 2300, S/D = 0.251). 
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Figure 2.21 . — Ensemble averaged velocity profile at r/D = 0.93 for Case 
lb (disc space flow, Va = 2300, S/D = 0.251). 


NAS A/CR— 2006-2 14131 


66 



Figure 2.22. — Ensemble averaged velocity profile at r/D = 1.05 for Case 
lb (disc space flow, Va = 2300, S/D = 0.251). 
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Figure 2.23. — Ensemble averaged velocity profile at r/D = 1.17 for Case 
lb (disc space flow, Va = 2300, S/D = 0.251). 
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Figure 2.24. — Ensemble averaged velocity profile at r/D = 1.30 for Case 
lb (disc space flow, Va = 2300, S/D = 0.251). 



Figure 2.25. — Turbulence intensity contour map at r/D = 0.68 for Case 
lb (disc space flow, Va = 2300, S/D = 0.251). 
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Figure 2.26. — Turbulence intensity contour map at r/D = 0.81 for Case 
lb (disc space flow, Va = 2300, S/D = 0.251). 
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Figure 2.27. — Turbulence intensity contour map at r/D = 0.93 for Case 
lb (disc space flow, Va = 2300, S/D = 0.251). 
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Figure 2.28. — Turbulence intensity contour map at r/D = 1 .05 for Case 
lb (disc space flow, Va = 2300, S/D = 0.251). 
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Figure 2.29. — Turbulence intensity contour map at r/D = 1 .17 for Case 
lb (disc space flow, Va = 2300, S/D = 0.251). 
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Figure 2.30. — Turbulence intensity contour map at r/D = 1 .30 for Case 
lb (disc space flow, Va = 2300, S/D = 0.251). 
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Figure 2.3 1 . — Ensemble averaged velocity profile for case Ha (tube flow, 
Va = 2300, SAD = 0.591). 
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Figure 2.32. — 77 contour map for case Ha (tube flow, Va = 2300, S/D = 
0.591). 
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Figure 2.33. — Flow visualization results in the degree range 0° < 0 < 15° for case Ha 
(tube flow, Va = 2300, S/D = 0.591). 
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Figure 2.34. — Flow visualization results in the degree range 15° < 0 < 60° for case Ha 
(tube flow, Va = 2300, S/D = 0.591). 
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Figure 2.35. — Flow visualization results in the degree range 60° < 0 < 90° for case Ha 
(tube flow, Va = 2300, S/D = 0.591). 
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Figure 2.36. — Flow visualization results in the degree range 90° < 0 < 150° for case Ha 
(tube flow, Va = 2300, S/D = 0.591). 
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Figure 2.37. — Flow visualization results in the degree range 180° < 0 < 240° for case 
Ha (tube flow, Va = 2300, S/D = 0.591). 
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Figure 2.38. — Flow visualization results in the degree range 240° < 9 < 330° for case 
Ha (tube flow, Va = 2300, S/D = 0.591). 



Figure 2.39. — The main features of the blowing stroke as revealed by flow 
visualization techniques for case lib (disc space flow, Va = 2300, S/D = 0.591). 
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Figure 2.40. — The main features of the suction stroke as revealed by flow visualization 
techniques for case lib (disc space flow, Va = 2300, S/D = 0.591). 



0.12 
0.1 
0.08 
3 " 0.06 

0.04 

0.02 

0 

1 


Figure 2.41 . — Velocity profile at r/D = 1 .30 for case lib (disc space flow, 
Va = 2300, S/D = 0.591). 
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Figure 2.42. — Turbulence intensity contour map at r/D = 1 .30 for case 
lib (disc space flow, Va = 2300, S/D = 0.591). 
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Figure 2.43 . — Velocity profile at r/D = 1 . 1 7 for case lib (disc space flow, 
Va = 2300, S/D = 0.591). 
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Figure 2.44. — Turbulence intensity contour map at r/D = 1 .17 for case 
lib (disc space flow, Va = 2300, S/D = 0.591). 
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Figure 2.45. — Velocity profile at r/D = 1 .05 for case lib (disc space flow, 
Va = 2300, S/D = 0.591). 
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Figure 2.46. — Turbulence intensity contour map at r/D = 1 .05 for case 
lib (disc space flow, Va = 2300, S/D = 0.591). 



Figure 2.47. — A sample of flow visualization data for case lib (disc space flow, Va = 
2300, S/D = 0.591). 
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Figure 2.48. — Velocity profile at r/D = 0.93 for case lib (disc space flow, 
Va = 2300, S/D = 0.591). 



Figure 2.49. — Turbulence intensity contour map at r/D = 0.93 for case 
lib (disc space flow, Va = 2300, S/D = 0.591). 
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Figure 2.50. — Velocity profile at r/D = 0.81 for case lib (disc space flow, 
Va = 2300, S/D = 0.591). 



Figure 2.5 1 . — Turbulence intensity contour map at r/D = 0.81 for case 
lib (disc space flow, Va = 2300, S/D = 0.591). 
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Figure 2.52. — Velocity profile at r/D = 0.68 for case lib (disc space flow, 
Va = 2300, S/D = 0.591). 



Figure 2.53. — Turbulence intensity contour map at r/D = 0.68 for case 
lib (disc space flow, Va = 2300, S/D = 0.591). 
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Figure 2.54. — A flow visualization result for 45 < 9 < 60 for case lib (disc 
space flow, Va = 2300, S/D = 0.591). 
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Figure 2.55. — Shortly after TDC for Case Ilia (tube flow, Va = 5400, S/D = 0.251). 

The leftward facing (long) arrows represent particles that are rapidly accelerated, the 
(shorter) curved arrows represent particles that maintained momentum in the opposite 
direction of the previous half-cycle and are still reversing. 
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Figure 2.56.-18° < 6 < 72° after TDC for Case Ilia (tube flow, Va = 5400, S/D = 
0.251). Growing separation bubbles with a low turbulence core. 
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Figure 2.57.-72° < 0 < 126° after TDC for Case Ilia (tube flow, Va = 5400, S/D = 
0.251) . Continued bubble growth with oscillating inlet conditions. 
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Figure 2.58.— 26° < 0 < BDC for Case Ilia (tube flow, Va = 5400, S/D = 0.251). Fully 
turbulent core flow, higher growth rate of ‘top’ separation bubble in cross-stream 
direction. Within the separation bubbles, random turbulent behavior is observed. 
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Figure 2.59. — Low-turbulence slug (left) pushing highly-turbulent slug out of the 
cylinder (right). This regime applies to all crank angles from BDC to TDC for Case 
Ilia (tube flow, Va = 5400, S/D = 0.25 1). 
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Figure 2.60. — End of cycle for Case Ilia (tube flow, Va = 5400, S/D = 0.251). Note that 
there is a non-uniform velocity distribution in the approximately stagnant fluid. The 
fluid in the upper portion of the cylinder maintains some residual momentum. It 
slowly displaces the fluid below. 
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Figure 2.61 . — Velocity profile for Case Ilia (tube flow, Va = 5400, 
S/D = 0.251). 
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Figure 2.62. — Turbulence intensity contour map for Case Ilia (tube flow, 
Va = 5400, S/D = 0.25 1). 
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Figure 2.63. — Flow pattern for small crank angles of the suction stroke 
for case Illb (disc space flow, Va = 5400, S/D = 0.251). 
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Figure 2.64. — Flow pattern for intermediate crank angles for Case Illb (disc space 
flow, Va = 5400, S/D = 0.251). 
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Figure 2.65. — Flow pattern for crank angles near maximum piston speed for Case Illb 
(disc space flow, Va = 5400, S/D = 0.251). 
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Figure 2.66. — Flow pattern during the blowing stroke period of strong 
acceleration (190° < 0 < 310°) for Case Illb (disc space flow, Va = 5400, 
S/D = 0.251). 
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Figure 2.67. — The blowing stroke for crank angles greater than 310° for Case Illb (disc 
space flow, Va = 5400, S/D = 0.251). 
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Figure 2.68. — Velocity profile at r/D = 0.68 for Case Illb (disc space 
flow, Va = 5400, S/D = 0.251). 



Figure 2.69. — 77 contour plot at r/D = 0.68 for Case Illb (disc space 
flow, Va = 5400, S/D = 0.251). 
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Figure 2.70. — Velocity profile at r/D = 0.81 for Case Illb (disc space 
flow, Va = 5400, S/D = 0.251). 



Figure 2.7 1 . — 77 contour plot at r/D = 0.8 1 for Case Illb (disc space 
flow, Va = 5400, S/D = 0.251). 
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Figure 2.72. — Velocity profile at r/D = 0.93 for Case Illb (disc space 
flow, Va = 5400, S/D = 0.251). 



Figure 2.73. — 77 contour plot at r/D = 0.93 for Case Illb (disc space 
flow, Va = 5400, S/D = 0.251). 
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Figure 2.74. — Velocity profile at r/D = 1 .05 for Case Illb (disc space 
flow, Va = 5400, S/D = 0.251). 



Figure 2.75. — 77 contour plot at r/D = 1.05 for Case Illb (disc space 
flow, Va = 5400, S/D = 0.251). 


NAS A/CR— 2006-2 14131 


93 



x/S e 


Figure 2.76. — Velocity profile at r/D =1.17 for Case Illb (disc space 
flow, Va = 5400, S/D = 0.251). 



Figure 2.77. — 77 contour plot at r/D = 1.17 for Case Illb (disc space 
flow, Va = 5400, S/D = 0.251). 
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Figure 2.78. — Velocity profile at r/D = 1 .30 for Case Illb (disc space 
flow, Va = 5400, S/D = 0.251). 



Figure 2.79. — 77 contour plot at r/D = 1 .30 for Case Illb (disc space 
flow, Va = 5400, S/D = 0.251). 
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Figure 2.80. — Wall Shear Stress and the associated error at r/D = 0.68 
for Case lib (Re max = 7600, Va = 2300, S = 127 mm, S/D = 0.591). 




Figure 2.81 . — Wall Shear Stress and the associated error at r/D = 0.81 
for Case lib (Re max = 7600, Va = 2300, 5 = 127 mm, S/D = 0.591). 
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Figure 2.82. — Wall Shear Stress and the associated error at r/D = 0.93 
for Case lib (Re max = 7600, Va = 2300, S = 127 mm, S/D = 0.591). 




Figure 2.83. — Wall Shear Stress and the associated error at r/D = 1.05 
for Case lib (Re max = 7600, Va = 2300, S = 127 mm, S/D = 0.591). 
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Figure 2.84. — Wall Shear Stress and the associated error at r/D =1.17 
for Case lib (Re max = 7600, Va = 2300, S = 127 mm, S/D = 0.591). 




Figure 2.85. — Wall Shear Stress and the associated error at r/D = 1.30 
for Case lib (Re max = 7600, Va = 2300, 5 = 127 mm, S/D = 0.591). 
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gure 2.86. — Wall Shear Stress and the associated error at r/D = 0.68 
for Case lb (Re max = 7600, Va = 2300, S = 54 mm, S/D = 0.25 1). 




Figure 2.87. — Wall Shear Stress and the associated error at r/D = 0.81 
for Case lb (Re max = 7600, Va = 2300, S = 54 mm, S/D = 0.25 1). 
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Figure 2.88. — Wall Shear Stress and the associated error at r/D = 0.93 
for Case lb (Re max = 7600, Va = 2300, 5 = 54 mm, S/D = 0.25 1). 




Figure 2.89. — Wall Shear Stress and the associated error at r/D = 1.05 
for Case lb (Re max = 7600, Va = 2300, S = 54 mm, S/D = 0.25 1). 
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Figure 2.90. — Wall Shear Stress and the associated error at r/D =1.17 
for Case lb (Re max = 7600, Va = 2300, 5 = 54 mm, S/D = 0.25 1). 



Figure 2.91 . — Wall Shear Stress and the associated error at r/D = 1 .30 
for Case lb (Re max = 7600, Va = 2300, 5 = 54 mm, S/D = 0.25 1). 
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gure 2.92. — Wall Shear Stress and the associated error at r/D = 0.68 
for Case Illb (Re max = 7600, Va = 2300, S = 54 mm, S/D = 0.251). 




Figure 2.93. — Wall Shear Stress and the associated error at r/D = 0.81 
for Case Illb (Re max = 7600, Va = 2300, 5 = 54 mm, S/D = 0.251). 


NAS A/CR— 2006-2 14131 


102 


RSS error T waii [ N/m 1 * RSS error T waii [ N/m 1 




gure 2.94. — Wall Shear Stress and the associated error at r/D = 0.93 
for Case Illb (Re max = 7600, Va = 2300, 5 = 54 mm, S/D = 0.251). 




Figure 2.95. — Wall Shear Stress and the associated error at r/D = 1.05 
for Case Illb (Re max = 7600, Va = 2300, 5 = 54 mm, S/D = 0.251). 
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Figure 2.96. — Wall Shear Stress and the associated error at r/D =1.17 
for Case Illb (Re max = 7600, Va = 2300, 5 = 54 mm, S/D = 0.251). 




Figure 2.97. — Wall Shear Stress and the associated error at r/D = 1.30 
for Case Illb (Re max = 7600, Va = 2300, 5 = 54 mm, S/D = 0.251). 
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Figure 2.98. — An approximate picture of the flow around the probe and probe support. 
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Figure 2.99. — Near wall measurements and wall conduction effects. 
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Figure 2.100. — An example of wall conduction effects. 



Figure 2.101 . — Flow leakage through access hole. 
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Figure 2. 102. — Velocity disturbance near x/S = 1 . 
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3. A Comparison Between Oscillatory Flow and Unidirectional Flow 

3.1 Introduction 

One may ask, given the data presented in sections 1 and 2, “How similar is oscillatory flow to steady 
unidirectional flow?” The impetus is that there is a large body of literature dedicated to impinging jet 
investigations, whereas, there is limited data available for oscillating flows in Stirling engines. Thus, to 
answer this question, we draw a comparison between steady and oscillatory flow. Presented below are 
comparisons between three different operating points. 

In theory, one way to predict whether unidirectional results can be applied to oscillatory flow is by 
looking at the Valensi number, 


Va = 


coD£ 

4v 


The Valensi number is a ratio of a recovery time (the effect of viscosity, i.e., if a fluid is quite viscous 
then it ‘recovers’ quickly - velocity gradients in the fluid will be quickly damped and the fluid will 
‘recover’) to that of the rotational period (oscillation period). For example, if the viscosity of an 
oscillating fluid is very small compared to the oscillation period, then the Valensi number is large. 
Physically speaking, velocity gradients within the fluid are not equalizing throughout the cycle - in other 

words, the fluid is not behaving like it is under steady state conditions. The quantity D h (the hydraulic 

diameter) appears in the Valensi number. This length scale tries to capture the length of a velocity 
gradient - that is, how far a disturbance must diffuse. 

The Valensi numbers presented in table 2.1 of section 2 are based on crankshaft rotational frequency, 
the viscosity of air and the tube diameter. This leads to a question: Is it valid to use tube-based Valensi 
numbers in the disc space? The answer is no. Another question arises in light of the very turbulent nature 
of the flow in portions of the test section during certain portions of the oscillation cycle: Is it appropriate 
to use the molecular viscosity of air to compute a Valensi Number? The answer to that is more 
complicated. Niu et al. 1 have shown that, in a regenerator matrix, turbulent eddy viscosity can be as much 
a 500 times as large as molecular viscosity. The implications of the work by Niu et al. 1 on the 
investigations under consideration are unclear - it is however something to contemplate. 

We present below a table (table 3.1) with Valensi numbers based on the length scale appropriate to 
the area of the test section under consideration. 


TABLE 3.1.— ADJUSTED VALENSI NUMBERS 


Case 

Location 

rpm 

S/D 

S 

(mm) 

Lengthscale 

Va 

Laminar 

Va 

Turbulent 

I.a 

Tube 

30 

0.251 

54 

D 

2300 

23 

I.b 

Channel 

30 

0.251 

54 

S 

580 

6 

Il.a 

Tube 

30 

0.591 

127 

D 

2300 

23 

Il.b 

Channel 

30 

0.591 

127 

S 

3213 

32 

III. a 

Tube 

70 

0.251 

54 

D 

5400 

54 

IILb 

Channel 

70 

0.251 

54 

S 

1361 

14 


Turbulent Valensi numbers were calculated by assuming that turbulent eddy viscosity is 100 times as 
large as molecular viscosity. This factor of 100 is strictly a guess at the effect of turbulent transport - 
further research is necessary to arrive at a more concrete value. 
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We see that all of the Valensi numbers in the disc space are much greater than unity regardless of the 
choice of viscosity . Thus, we expect that for the cases presented below, there will be little to no 
correlation between oscillatory and unidirectional results. 

3.2 Nomenclature 


Symbol 

Description 

Units 

D 

Delivery tube diameter 

[m] 

Re 

Reynolds Number 

[dimensionless] 

r 

Radial Variable 

[m] 

r/D 

Non-Dimensional radial variable 

[dimensionless] 

rpm 

revolutions per minute 

[1 /min] 

S 

Disc Spacing 

[m] 

Va 

Valensi number 

[dimensionless] 

X 

Axial variable 

[m] 

x/S 

Non-Dimensional axial variable 

[dimensionless] 

0 

Crank Angle 

[degrees] 

CO 

Rotational speed 

[1/min] 


3.3 Re = 5989, S = 127 mm, co = 30 rpm, Va - Laminar = 3213, Va - Turbulent = 32 

In this case we look at the disc space flow and draw a comparison between the unidirectional case 
Re = 7600, S = 127 mm (Case 1 in section 1), and the oscillatory flow case, co = 30 rpm, S = 127 mm 
(Case lib in section 2). We begin with a note; in the unidirectional investigation, Case 1 is affiliated with 
a Reynolds number of 7600. This is not exactly correct (as explained in section 1 subheading 1 .7. 1 ). In 
section 1, heading 1.7.1, Reynolds numbers were based on tube centerline velocities. At present, this is 
not the correct course of action; we must observe a mass balance. To this end, we base Reynolds number 
on the bulk average velocity across the delivery tube under steady flow conditions. If the Reynolds 
number is based on bulk average velocity, then under unidirectional conditions, Re = 5989. 

Now, under oscillatory flow, recall that in Case lib of section 2, the cited maximum Reynolds 
Re max number is 7600 (see table 2.1, section 2). Recall also that Reynolds number in the oscillatory 
studies is based on piston velocity, and, hence, is driven by a sinusoidal pressure gradient. Thus, to match 
a unidirectional Reynolds number of 5989, two crank angles are possible, one while the flow is 
accelerating on the blowing stroke, and one while the flow is decelerating on the blowing stroke. These 
two angles are denoted as 0j = 232° and 0 9 = 308° . Note that the maximum piston velocity occurs at 
0 = 270. The results of this comparison are presented in figures 3.1 and 3.2. 

Figure 3.1 shows that there is little to no correlation between oscillatory and unidirectional flow 
results in terms of both profile shape and velocity magnitude. The lack of agreement clearly shows the 
influence of unsteady effects in such flows. As in the case of the accelerating flow (fig. 3.1), figure 3.2 
shows no correlation with decelerating flow. Hence, figure 3.2 provides additional verification that 
unsteady flow effects govern the flow field. This is not surprising in light of the Valensi number 
associated with this case. 

Figure 3.3 shows the effect of flow acceleration on the average velocity profiles. Both angles 
correspond to a Reynolds number Re = 5989. However, at 0, , the piston acceleration is 62% of the peak 
value, whereas, at 0 2 , the piston acceleration is -62% of it’s peak value (i.e., the piston is decelerating). 
The effect that piston acceleration has on fluid velocity is remarkable. 


We see that one of the Valensi numbers based on turbulent viscosity is 6.0. We may hope that this case is almost 
quasi steady. We shall see. 
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The results of figures 3.1, 3.2, and 3.3 lead to a comment; it is now apparent that although the piston 
follows an exact sinusoid, there is no reason to expect the same from any given point in the flow field. 
This is clear to one who looks at the flow visualization results from section 2 (figs. 3.4 and 3.5). 

Figures 3.4 and 3.5 show the drastic change in the flow field across the maximum piston velocity 
( 0 = 270). The flow velocity definitely does not follow the piston velocity. 

3.4 Re = 14937, S = 54 mm, co = 70 rpm, Va - Laminar = 1361, Va - Turbulent = 14 

Here we compare unidirectional results from Case 4 of section 1 (RE = 17,700 S = 54 mm) to 
oscillatory results from Case Illb of section 2 (to = 70 rpm, S = 54 mm). Once again, with regard to 
unidirectional flow, we have a Reynolds number based on average velocity across the tube and tube 
diameter of Re = 14937, which corresponds to two crank angles under oscillatory flow. The crank angles 
of interest are 9j = 237.5°and 0-> = 302.5°. Presented in figure 3.6 is a comparison between unidirectional 
results and oscillatory results at 9j = 237.5°. We see that the profiles have similar shape but the velocity 
magnitudes are quite different. The wall jet in the unidirectional case is thinner and stronger. Turning to 
the second matching crank angle, we see (fig. 3.7) that there is general profile shape agreement and the 
velocity magnitudes seem to compare quite favorably in the wall jet region. For the sake of comparison, 
the two crank angles are plotted against each other in figure 3.8. We see that for 0 7 = 302.5°, the wall jet 
is stronger and the flow may be starting to reverse near the back wall. It seems that the fluid gains 
momentum as the blowing stroke proceeds, hence, stronger jets are seen at larger crank angles. The 
importance of the Valensi number is underscored in this case; figure 3.7 may suggest that in this case 
unidirectional results may be used to predict oscillatory behavior. However, the large Va value associated 
with this case, and the information presented in figure 3.8 suggest that this is not the truth. 

3.5 Re = 5839, S = 54 mm, © = 30 rpm, Va - Laminar = 580, Va - Turbulent = 6 

This is a comparison between Case 2 of section 1 (RE = 7600, S = 54 mm) and case lb of section 2 
(co = 30 rpm, S = 54 mm). With a unidirectional bulk average Reynolds number of 5839, the matching 
crank angles under oscillatory flow are 9j = 230° and 9? = 310°. Presented in figure 3.9 is a comparison 
between unidirectional results and 0j = 230° . Once again, the overall correlation between the two flow 
regimes is poor. The profile for the oscillatory case at r/D of 0.930 is remarkable; in this case the wall jet 
dies very rapidly in the radial direction. Looking at the other crank angle of interest (fig. 3.10) we see that 
for the smaller of the two radial locations, under oscillatory flow, the profile in the jet region surpasses 
that of the unidirectional flow profiles. In contrast, for the larger radial oscillatory flow profile, the 
velocity magnitude is considerably less than that of steady flow. Presented in figure 3.11 is a comparison 
between the two oscillatory cases. The effect of the addition of momentum and the utility of the Valensi 
number is reiterated. 


3.6 Conclusions 

The main and important result of this section is that Stirling engine designers should not try to apply 
steady-unidirectional impinging jet data to oscillatory flows for regions of the engine where the Valensi 
number is greater than unity. If an appropriate Va value were chosen (based on laminar diffusion or 
turbulent transport) then one would suspect that steady flow results would be applicable when the Valensi 
numbers are sub-unity. The problem lies in the determination of the Valensi number. In the test section 
under consideration, under oscillatory flow, highly turbulent regions were generated and stimulated by the 
oscillatory nature of the flow (For example, the separation zone that forms on the back wall and its re- 
ingestion into the tube space on the next half cycle). Thus, to apply unidirectional results, information 
would be required about the turbulent transport within the test section or engine for the purpose of 
determining turbulent viscosity. Further investigation into this topic is beyond the scope of this work. 
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3.7 Figures 



x/S 


Figure 3.1 . — A comparison of unidirectional results (Re = 5989) with oscillatory work 
at 0j = 232° => Re = 5989 . 
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0.35 



Figure 3.2. — A Comparison of unidirectional results (Re = 5989) with oscillatory work 
at 0 2 = 308° => Re = 5989 . 
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Figure 3.3. — A comparison between 0j = 232°, 0 2 = 308°. 
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Figure 3.4. — Flow Visualization for 240 < 0 < 255 (same piston velocity as that of fig. 3.5). 
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Figure 3.5. — Flow visualization for 285 < 0< 300(same piston velocity as that of fig. 3.4). 
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r/D = 0.827 (uni-d) 
r/D = 0.945 (uni-d) 
r/D = 0.807 (237.5 osc) 
r/D = 0.930 (237.5 osc) 



x/S 


Figure 3.6. — A comparison between steady unidirectional flow and oscillatory flow at 
0| = 237.5° with matching Reynolds numbers. 
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Figure 3.7. — A comparison between steady unidirectional flow and oscillatory flow at 
09 = 302.5° with matching Reynolds numbers. 
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r/D = 0.807 (237.5) 
r/D = 0.930 (237.5) 
r/D = 0.807 (302.5) 
r/D = 0.930 (302.5) 



0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 


x/S 


Figure 3.8. — A comparison between two crank angles ( 0 ( = 237.5°and 0 2 = 302.5°) 
with matching Reynolds numbers on the blowing stroke. 
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Figure 3.9. — A comparison between steady unidirectional flow and oscillatory flow at 
0| = 230° with matching Reynolds numbers. 
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Figure 3.10. — A comparison between steady unidirectional flow and oscillatory flow at 
09 = 310° with matching Reynolds numbers. 
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Figure 3.11 . — A comparison between two crank angles ( 0, = 248° and 0 2 = 292°) with 
matching Reynolds numbers on the blowing stroke. 
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4. Comparisons Between Experimentally 
Gathered Results and CFD Generated Data 


4.1 Introduction 

As has been alluded to several times prior to this point, one of the main reasons that this project was 
undertaken was code validation support. In this section, we make comparisons between the data presented 
in sections 1 and 2 and numerically simulated CFD results. We state now that the CFD results that will be 
presented below were not generated at the University of Minnesota. The computational results were 
generated at Cleveland State University. Before we proceed we wish to give credit to the following 
people at Cleveland State: Mounir B. Ibrahim, Principal Investigator, Zhiguo Zhang, Graduate Research 
Assistant and Sundeep Kembhavi, Graduate Research Assistant. 

4.2 Nomenclature 


Symbol 

Description 

Units 

Re 

Reynolds Number 

[dimensionless] 

S 

Disc Spacing 

[m] 

u 

x-direction velocity component 

[m/s] 

II u II 

/ 2 2 

Velocity magnitude ( II u 11= W + v ) 

[m/s] 

Va 

Valensi Number 

[dimensionless] 

0 

Crank Angle 

[degrees] 

CO 

Rotational Frequency 

[1 /min] 


4.3 Oscillatory Investigations 


4.3.1 A Comparison of Experimental Flow Visualization Results to CFD Generated Data for the 
Case Re max = 7600, S = 127 mm, Va = 2300, co = 30 rpm 

A comparison of the flow field in the disc space was made for an entire cycle. Velocity vectors and 
streamlines are plotted in the computational results (see fig. 4.1). Figure 4.1 also shows the area where the 
comparison with flow visualization is made. The flow visualization technique was explained in section 2, 
section 2.5.1. We will not repeat that discussion but will comment that the flow visualization pictures are 
ensemble-averaged results. The comparison is presented in figures 4.2 through 4.25. Note that the 
comparisons are made every 15° of crank angle. In general, the code seems to very accurately capture 
both the sink flow from the back wall on the suction stroke (0° < 0 < 180°) and the wall jet on the 
blowing stroke 180° < 0 < 360°. However, the recirculation zones are not described perfectly. For 
example, the code identifies eddies present in the flow. However, the locations of the eddies are not in 
accordance with the flow visualization results. In other words the computation is not accurately capturing 
the recirculation features of the flow. This is not surprising in the least when one considers that the 
computational pictures shown in this section were developed using no turbulence model. In other words, 
the computations were run without turbulence effects. 
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4.3.2 A Comparison of Experimentally Derived Velocity Data to CFD Generated Values for the 
Case Re max = 7600, S = 127 mm, Va = 2300, to = 30 rpm 

Presented in figures 4.26 through 4.33 are a comparison of CFD generated results vs. experimental 
data for oscillatory flow. Figures 4.26 through 4.29, focus on the crank angle 0 = 246.5° and each figure 
steps out (in order of appearance) in radial dimension, r/D = 0.473, 0.591, 0.827, and 0.945. The reader 
may notice that these radial locations are not the same as those discussed in section 2. The explanation is 
that the data for this comparison were generated specifically for two papers by Adolfson et al., 1 and 
Ibrahim et al. 2 Flowever, outside of the differences in radial locations, the methods for data collection and 
processing were for all practical purposes identical to those discussed in section 2. Figure 4.26 presents a 
comparison between experimental data and CFD generated results using two turbulence closure models; a 
laminar model and a K-to model. It seems that for this radial location (r/D = 0.473) the K-oo model does a 
better job of capturing the flow features near the target wall (small x/S). One can see that for x/S > 0.3, the 
models start to diverge from the experimental data. This is not entirely surprising in the sense that for 
larger x/S, the probe was (and hence the experimental data were) in a highly turbulent zone and the hot- 
wire measurements were not accurate. Thus, one should not expect a match in this area (for 77 plots for 
this case, see Adolfson et al. 1 ). Moving radially outward to r/D = 0.591, we see (fig. 4.27) that the laminar 
model does a better job of approximating the flow for 0 < x/S < 0.3. Again, for x/S > 0.3, turbulence 
quantities are too high to expect good hot wire measurements. Figures 4.28 and 4.29 present comparisons 
at r/D = 0.827 and 0.945, respectively. There is little correlation between either of the models and the 
experimental data at these radial locations. 

Figures 4.30 through 4.34 present results at a different crank angle, namely, 0 = 293.5°. The figures 
present results at the radial locations r/D = 0.473, 0.591, 0.827, and 0.945. Figures 4.30 and 4.31 show 
that the K-oo model accurately captures the velocity profile shapes but not the velocity magnitudes at the 
radial locations r/D = 0.473 and r/D = 0.591. Moving radially outward, figures 4.32 and 4.33 present 
comparisons between experiment and computation at the radial locations r/D = 0.827 and r/D = 0.945, 
respectively. At these radial locations, both models seems to capture profile shape reasonably well, but 
neither captures velocity magnitude. 

One further thing to notice about the CFD results is that they were presented at the crank angles 0 = 
293.5° and 0 = 246.5°. This corresponds to the cases discussed in section 3 where the accelerating and 
decelerating flow cases for the same Reynolds number are discussed. The computational model captures 
the effects discussed in section 3 (fig. 4.26 vs. 4.30, fig. 4.27 vs. 4.3 1, fig. 4.28 vs. 4.32, and fig. 4.29 vs. 
4.33). 


4.4 Unidirectional Investigations 

4.4.1 A Comparison of Experimental Flow Visualization Results to CFD Generated Data for the 
Case Re = 7600, S = 127 mm 

In section 1, subheading 1.6.1 we presented steady unidirectional flow visualization results. We now 
compare these experimental flow visualization results to CFD generated streamlines. For the readers’ 
information, this was originally presented by Ibrahim et al. 2 In figure 4.34 two models are presented; the 
leftmost figure plots streamlines as generated by a laminar flow model and the middle figure presents 
streamlines as generated with the same flow equation solver but using a K-co turbulence closure model. 
Quoting Ibrahim et al. 2 “As the flow turns into the channel [disc space], a large-scale vortex is created 
between the two discs. The flow accelerates from the stagnation zone, with additional acceleration due to 
the vortex generated between the two discs. Near r/R = 1 [r/D = 0.5] (r is the radial location from the axis 
of symmetry and R is the tube radius) the flow starts to decelerate due to the increase in flow area. The 
CFD results using the K-oo model seem to predict these flow features very well. On the other hand, the 
laminar flow model shows only a thin flow layer near the right [target] disc and not enough mixing to 
spread the flow across the channel.” 
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4.4.2 A Comparison of Experimentally Derived Velocity Data to CFD Generated Values for the 

Case Re = 7600, S = 127 mm 

Here, we present a comparison between unidirectional, steady flow anemometry measurements and 
CFD generated predictions using laminar and turbulent K-to closure models. In figure 4.35 results are 
presented at the radial location r/D = 0.473. One thing to note is that a CFD program will produce 
velocity components u=< u,v> , whereas the hotwire only measures velocity normal to the wire, 

II u ll=< u,v > ° < u,v >= (u + v ) , in this case. Consequently, the magnitudes of the CFD generated 

velocity components had been computed and plotted in figure 4.35. Looking at figure 4.35, one sees that 
in the near-wall region (small x/S), the K-co model compares well with experimental data. It is also 
apparent that for larger x/S, the laminar model is a closer approximation. Moving on to figure 4.36, a 
comparison is made at the radial location r/D = 0.591. It is apparent that for x/S < 0.5, neither model 
accurately captures velocity magnitude ||u||. However, the laminar model quite accurately captures the 
profile shape. Regarding figures 4.37 and 4.38, the plots are somewhat different. In these two plots, the 
radial velocity component u, as computed by the CFD routine, is plotted and compared to the hot wire 
measured speed, ||m||. The rational for this is that for the larger radial locations (r/D = 0.827 and r/D = 

0. 945, presented in figs. 4.37 and 4.38, respectively) the flow is mainly in the radial direction. This idea is 
supported by figure 1 .9 of section 1 . For larger radial locations, the flow is in one of two regimes; either it 
is in a wall jet that is moving radially outward, or it is in a separation bubble where hot wire derived 
measurements are not reliable. Regarding figure 4.37, it can be seen that the K-co model does a good job 
of approximating the flow near x/S = 0. Once again, the laminar model does a decent job of simulating the 
profile shape, but misses in magnitude. Moving radially out to r/D = 0.941 (fig. 4.38), the same type of 
behavior is observed. 

With regards to figures 4.35 through 4.38 we quote Ibrahim et al. 2 “Figure 3a [fig. 4.35] is for a radial 
location of r = 0.1016 m [r/D = 0.473].. .Here, near the wall region, the K-co model is in closer agreement 
with the experiment. The laminar flow is in a better agreement with the data a distance away from the 
wall. Figure 3b [fig. 4.36] shows a similar plot to the one shown in 3a but at r = 0.127m [r/D = 0.591]. 

The laminar flow seems to be in overall better agreement with the data than the K-co model results. 
Figures 3c and d [figs. 4.37 and 4.38] show the CFD data for the radial velocity component v [u]. Here 
we see the predictions using the K-co model are in closer agreement with the experiment. This is 
consistent with the discussion shown earlier for figure 2 [see the quote by Ibrahim et al. in section 4.5.1, 

1. e. the CFD generated streamlines vs. experimental flow viz.], where the vortex was predicted well using 
the K-co model and more spread of the velocity fields is noticed at higher radial locations.” 

4.5 Conclusions 

By and large, the CFD models do a very good job in approximating the bulk flow behavior within the 
test section under oscillatory and unidirectional flow conditions. The testament to this statement is the 
comparison with the flow visualization. However, in terms of accurately predicting velocities within the 
test section, the CFD routines fall short. It has been observed that under certain conditions, the laminar 
closure does a better job of capturing velocities, whereas, under different conditions, the K-co model 
yields superior results. The moral is that a CFD routine needs to be developed that has an adaptive closure 
model with regard to laminar or turbulent flow. Another point to note is that neither model captures 
transition to turbulence. This is problematic in the sense that transition influences the entire flowfield (by 
influencing the location of the separation bubble on the back disc). Thus, until transition is either 
accurately captured via the CFD routine (or accurately modeled for), one cannot expect accurate 
comparisons between CFD results and experimental data. 
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4.6 Figures 


In the figures that follow, you will see enlargements of CFD pictures 
generated by CSU. The enlargements were created from pictures 
like the one below. 



Figure 4.1 . — A sample of the type of data that can be generated with the CFD-ACE 
program used at Cleveland State University. Notice that streamlines and velocity 
vectors are plotted. 




Centerline of Cylinder 


Figure 4.2. — CFD generated results vs. flow visualization data for the oscillatory case 
Re ma x = 7600, Va = 2300, S = 127 mm during the crank angle range 0° < 0 < 15°. 
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Centerline of Cylinder 


Figure 4.3. — CFD generated results vs. flow visualization data for the 
oscillatory case Re max = 7600, Va = 2300, S = 127 mm during the crank angle 
range 15° < 0 < 30°. 




Centerline of Cylinder 


Figure 4.4. — CFD generated results vs. flow visualization data for the 
oscillatory case Re max = 7600, Va = 2300, 5 = 127 mm during the crank angle 
range 30° < 0 < 45°. 
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Centerline of Cylinder 


Figure 4.5. — CFD generated results vs. flow visualization data for the 
oscillatory case Re max = 7600, Va = 2300, S = 127 mm during the crank angle 
range 45° < 0 < 60°. 




Centerline of Cylinder 


Figure 4.6. — CFD generated results vs. flow visualization data for the 
oscillatory case Re max = 7600, Va = 2300, S = 127 mm during the crank angle 
range 60° < 0 < 75°. 
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Centerline of Cylinder 


Figure 4.7. — CFD generated results vs. flow visualization data for the 
oscillatory case Re max = 7600, Va = 2300, S = 127 mm during the crank angle 
range 75° < 0 < 90°. 




Centerline of Cylinder 


Figure 4.8. — CFD generated results vs. flow visualization data for the 
oscillatory case Re max = 7600, Va = 2300, S = 127 mm during the crank angle 
range 90° < 0 < 105°. 
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Centerline of Cylinder 


Figure 4.9. — CFD generated results vs. flow visualization data for the 
oscillatory case Re max = 7600, Va = 2300, S = 127 mm during the crank angle 
range 105° < 0 < 120°. 




Centerline of Cylinder 


Figure 4.10. — CFD generated results vs. flow visualization data for the 
oscillatory case Re max = 7600, Va = 2300, S = 127 mm during the crank angle 
range 120° < 0 < 135°. 
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Centerline of Cylinder 


Figure 4.1 1 . — CFD generated results vs. flow visualization data for the 
oscillatory case Re max = 7600, Va = 2300, S = 127 mm during the crank angle 
range 135° < 0 < 150°. 




Centerline of Cylinder 


Figure 4.12. — CFD generated results vs. flow visualization data for the 
oscillatory case Re max = 7600, Va = 2300, S = 127 mm during the crank angle 
range 150° < 0 < 165°. 
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Centerline of Cylinder 


Figure 4.13. — CFD generated results vs. flow visualization data for the 
oscillatory case Re max = 7600, Va = 2300, S = 127 mm during the crank angle 
range 165° < 0 < 180°. 




Centerline of Cylinder 


Figure 4.14. — CFD generated results vs. flow visualization data for the 
oscillatory case Re max = 7600, Va = 2300, S = 127 mm during the crank angle 
range 180° < 0 < 195°. 
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Centerline of Cylinder 


Figure 4.15. — CFD generated results vs. flow visualization data for the 
oscillatory case Re max = 7600, Va = 2300, S = 127 mm during the crank angle 
range 195° < 0 < 210°. 



Figure 4.16. — CFD generated results vs. flow visualization data for the 
oscillatory case Re max = 7600, Va = 2300, S = 127 mm during the crank angle 
range 210° < 0 < 225°. 
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Figure 4.17. — CFD generated results vs. flow visualization data for the 
oscillatory case Re max = 7600, Va = 2300, S = 127 mm during the crank angle 
range 225° < 0 < 240°. 
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Figure 4.18. — CFD generated results vs. flow visualization data for the 
oscillatory case Re max = 7600, Va = 2300, S = 127 mm during the crank angle 
range 240° < 0 < 255°. 
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Figure 4.19. — CFD generated results vs. flow visualization data for the 
oscillatory case Re max = 7600, Va = 2300, S = 127 mm during the crank angle 
range 255° < 0 < 270°. 



Figure 4.20. — CFD generated results vs. flow visualization data for the 
oscillatory case Re max = 7600, Va = 2300, S = 127 mm during the crank angle 
range 270° < 0 < 285°. 
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Figure 4.21 . — CFD generated results vs. flow visualization data for the 
oscillatory case Re max = 7600, Va = 2300, S = 127 mm during the crank angle 
range 285° < 0 < 300°. 
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Centerline of Cylinder 

Figure 4.22. — CFD generated results vs. flow visualization data for the 
oscillatory case Re max = 7600, Va = 2300, S = 127 mm during the crank angle 
range 300° < 0 < 315°. 
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Centerline of Cylinder 


Figure 4.23. — CFD generated results vs. flow visualization data for the 
oscillatory case Re max = 7600, Va = 2300, S = 127 mm during the crank angle 
range 315° < 0 < 330°. 



Figure 4.24. — CFD generated results vs. flow visualization data for the 
oscillatory case Re max = 7600, Va = 2300, S = 127 mm during the crank angle 
range 330° < 0 < 345°. 
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Figure 4.25. — CFD generated results vs. flow visualization data for the 
oscillatory case Re max = 7600, Va = 2300, S = 127 mm during the crank angle 
range 345° < 0 < 360°. 



Figure 4.26. — CFD vs. experimental data under oscillatory flow when 0 = 246.5° at the 
radial location r/D = 0.473. 
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Figure 4.27. — CFD vs. experimental data under oscillatory flow when 0 = 246.5° at the 
radial location r!D = 0.591 . 



Figure 4.28. — CFD vs. experimental data under oscillatory flow when 0 = 246.5° at the 
radial location r/D = 0.827. 
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Figure 4.29. — CFD vs. experimental data under oscillatory flow when 0 = 246.5° at the 
radial location r/D = 0.945. 
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Figure 4.30. — CFD vs. experimental data under oscillatory flow when 0 = 293.5° at the 
radial location r/D = 0.473. 
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Figure 4.3 1 . — CFD vs. experimental data under oscillatory flow when 0 = 293 .5° at the 
radial location r!D = 0.591 . 



Figure 4.32. — CFD vs. experimental data under oscillatory flow when 0 = 293.5° at the 
radial location r/D = 0.827. 
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Figure 4.33. — CFD vs. experimental data under oscillatory flow when 0 = 293.5° at the 
radial location r!D = 0.945. 




Figure 4.34. — A comparison between CFD generated streamlines and experimental 
flow visualization under unidirectional flow when Re = 7600, S = 127 mm. 
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Figure 4.35. — CFD vs. experimental data for the case Re = 7600, S = 127 mm at the 
radial location r/D = 0.473. 



x/S 

Figure 4.36. — CFD vs. experimental data for the case Re = 7600, S = 127 mm at the 
radial location r/D = 0.591. 
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Figure 4.37. — CFD vs. experimental data for the case Re = 7600, S = 127 mm at the 
radial location r!D = 0.827. 
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Figure 4.38. — CFD vs. experimental data for the case Re = 7600, S = 127 mm at the 
radial location r/D = 0.945. 
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Appendix A — Hot Wire Anemometry Basics 

A.l Introduction 

In this appendix, we attempt to familiarize the reader with a hot wire anemometer. To this end, we 
provide a picture of an anemometry system, list the quantities it measures, and explain how it works. 
Starting out simply, the hot wire anemometer is a device that can be used to measure the product of 
density and velocity (p«). From this quantity, under certain assumptions (i.e., constant, or 
deterministically variable density and velocity) the velocity, temperature, or density of a flowing fluid can 
be measured. 

Before we proceed, we present a picture of a hot wire anemometer system in figure A. 1 . Shown in the 
figure is an anemometer, a probe cable, a probe support and a probe. The general idea is that the user 
exposes the probe to a flow. The probe then measures the velocity of the fluid around it. The actual 
mechanics of the measurement will be discussed shortly, but before that, we introduce one more figure, 
figure A.2. This figure shows a typical hot-wire probe. The main features of the probe are a wire, support 
prongs, a probe body, and leads. The main object of concern is the wire. It is generally made of platinum, 
tungsten or a platinum-tungsten alloy. The diameter of the wire is generally on the order of 5 pm. 


A.2 Nomenclature 


Symbol 

Description 

Units 

D w 

Wire diameter 

M 

e 

Voltage differential across a Wheatstone bridge 

[ R ] 

h 

Average convective heat transfer coefficient 

Wlm 2 °K 

i 

Electrical current 

[A] 

k 

Fluid thermal conductivity 

[W/moK] 

L 

Wire length 

M 

Nu 

Average Nusselt number 

[] 

Pr 

Prandtl number 

[] 

q' 

Heat transfer per unit length 

[W / m\ 

R 

Electrical resistance 

[«] 

Re w 

Wire diameter Reynolds number 

[] 

T 

1 amb 

Ambient temperature 


T 

1 w 

Wire temperature 


V 

Measured voltage 

[v] 

X 

Velocity magnitude (speed) 

[m / 5 ] 

r\ 

Variable resistance 

m 

CO 

Total probe resistance 

[n] 


A.3 Basic Anemometer Circuitry 

For the present discussion, the main thing to keep in mind is that the wire is a temperature sensitive 
resistive element. In other words, given a wire temperature, there is an associated resistance and vice 
versa. This concept will be of use shortly in the discussion of operating resistance; namely, in specifying 
an operating resistance, an operating temperature is implied. 
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Shown in figure A.3 is a Wheatstone bridge. The Wheatstone bridge is the fundamental circuit in an 
anemometer. This circuit is said to be in balance if there is a net flow of current and: 

e = 0 V. 

For the above equation to be satisfied, it must be the case that (referring to fig. A1.3): 

R R 

co T| 

In the case of a hot wire anemometer, we let r| be a variable resistor and let co be the “operating 
resistance” of the probe. “The probe [operating] resistance is a composite of the resistance of the sensor 
[at the operating temperature] plus the resistance of the support needles, probe body, leads and probe 
cable” (Lomas, 3 p. 85). 

If the probe resistance and the operating resistance are equal, then the bridge is in balance. Exposing a 
probe to a flow will cause a change in heat transfer via convection (to be explained shortly). As a result, 
its temperature and, hence, its resistance will change. The net result is that: 

e*0 V. 

This voltage information is supplied to a control system (a series of amplifiers that are not shown in 
the figure) which varies the current through the probe. The change in current induces a change in resistive 
dissipation. This raises or lowers the heating rate of the probe and returns it to the specified operating 
resistance; i.e., the bridge is in balance again. 

A.4 Physical Principles of an Anemometer Probe 

A thermal anemometer measures the cooling rate of a small cylinder in cross flow. From the 
measured cooling rate, a speed associated with the flowing fluid can be deduced. We set forth to explain 
the physics behind an anemometer probe in a somewhat heuristic manner. 1 

As a simple model, imagine an infinitely long cylinder. Imagine further that there is a one 
dimensional, constant property, flowing fluid with its velocity vector normal to the cylinder. We allow the 
fluid velocity to vary. Now, suppose that we can somehow heat the cylinder such that its temperature is 
uniform in space and time. Given these assumptions, we begin a “simple” hot-wire analysis. For clarity, 
the situation is presented pictorially in figure A.4. 

We begin by focusing on how thermal energy is lost from the cylinder. The most fundamental starting 
point is Newton’s law of cooling: 


q'= nD w h(T w - T amh ) . 

We introduce the fluid thermal conductivity and rearrange the above equation to yield: 

<?' hD » 


kn(T w - T amh ) 


- = Nu 


1 When we say heuristic, we are seeking to give the reader a quick and dirty derivation of King’s law that is easy to 
grasp. We have not followed “conventional” derivations to the letter of the law; instead we have made a tractable 
physical argument. For a more conventional derivation (which is beyond the scope of this work) the reader is 
referred to Lomas 3 (pp. 55-72). 
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Nusselt numbers are tabulated in a wide variety of basic texts on heat transfer. To find a correlating 
function, we cite Incropera and De Witt. 1 (p. 344): 

c i = Nu = C Re'” Pr 1/3 . 

kn(T w -T amh ) 

where C and m are given for various Re w ranges. 

Solving for the heat transfer from the cylinder per unit length, we have: 

q'= Ckn(T w - T amb ) Re” Pr 1/3 . 

We now turn our attention to the heating of the cylinder. In the case of a hot-wire anemometer, the 
cylinder (wire) is heated electrically. Recall that the resistance of the wire is a function of temperature; 
thus, if temperature is constant (which we demanded initially) then the resistance of the cylinder is 
constant. We now state Ohm’s law: 


V = iR, 

from which we write the heat generation per unit length as: 


<7 = 


XI 

RL 


We have the necessary components to relate heat transfer by convection to heat generated by resistive 
dissipation, viz. 


y2 

- = CkK(T w -T amb ) Re^Pr 1/3 . 

Grouping all constants and rewriting we have: 

V 2 = Cx m . 

Where C and m are functions of the Reynolds and Prandtl numbers. 

The above relation suggests that at zero velocity, the voltage output from the hot wire should be 
identically zero. This is invariantly not the case. At zero fluid velocity, there is always a non-zero 
anemometer voltage output. To explain this phenomenon, we turn to physical arguments. 

The basic concept is simple, if we were to have zero voltage at zero flow, then there would be no loss 
of thermal energy from the wire. Since the wire is less than perfectly isolated from its surroundings, there 
is a loss of energy. First of all, the hot wire is not of infinite extent; at each end there is a support prong. 
At these prongs, there will be a conduction loss (i.e., heat will flow from the wire to the prongs). Also, it 
is the case that most fluids become less dense as their temperature increases (water near the freezing 
point, being and exception). In any case, the wire is hot and the fluid around it is cool. The wire heats the 
fluid around it (via convection) and the fluid rises away from the wire. The effect is that the wire is 
continuously losing thermal energy via ‘natural convection.’ A third form of energy loss is due to radiant 
effects. A typical hot wire operating temperature is 250 °C. The effects of radiation at these low 
temperatures may be small, but not zero. 

The point is that at zero flow there is still heat transfer, and, hence, a current or voltage is needed to 
balance the Wheatstone bridge. This suggests that we may improve our above-derived result by adding a 
correction to it, namely: 


NAS A/CR— 2006-2 14131 


149 



V 2 = A + Cx n 

where A, C, and m are now unknown and considered to be constant. Thus, forms of heat transport other 
than forced convection are allowed and assumed to have an effect that is independent of velocity. Observe 
that at very low velocities, the main contribution to the voltage signal will come from these alternative 
modes. 

Replacing m with 0.5 yields a well-known result called King’s law (King 2 ). If m is taken as a given 
constant (not necessarily 0.5), the above form is very useful. It allows one to (theoretically) take two 
voltage signals, calibrate the anemometer to two corresponding velocity points and with no further data 
required, plot a slope-intercept line that describes a wide range of voltage-velocity points. When we say 
theoretically, we suggest that one must know the value of m to carry out a 2-point calibration. 

So, given a data set, how do we determine the value of ml Generally speaking, given a velocity- 
voltage data set, the procedure is to guess a value of m, plot V 2 vs. x n , fit a line to the result and check 
how well the linear fit approximates the data points. This process is carried out iteratively until a best fit is 
found. It turns out that the optimum value of m is 0.435 for the probes of the present study. As a benefit, 
when the best linear fit is found, both A and C are determined as well - this whole process is called 
calibration, which is the subject of appendix B. 


A.5 Figures 



Figure A.l. — An example of a hot-wire anemometry system. 
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Figure A.2. — A TSI model 1210 hot-wire probe. 



Figure A.3. — A schematic of a Wheatstone bridge. 
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Figure A.4. — A schematic of the assumptions necessary for a simple hot wire analysis. 
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Appendix B — Qualification, Error and 
Calibration of a Hot-Wire in a Low Speed Flow 

B.l Introduction 

Of particular interest in the investigations above is the behavior of a hot wire under low velocity flow 
conditions. In an oscillating flow, there are portions of the cycle where the flow has zero velocity. 
Furthermore, all of the velocities measured above, both under oscillatory and unidirectional flow 
conditions, were lower than 2 m/s (sometimes by an order of magnitude or more). The question of interest 
is when does the anemometer output drift from linear behavior (see appendix A if clarification of this 
statement is required). 


Symbol 

C ch 
Cl 

D w 

Dch 

DelP 

e 


8 

Gr w 

h 


i 

k 

L 

Nu 


[ atm 


P c h 

Pr 

q' 

R 


Rci w 

Re c/! 

Re w 


RTDR 


L amb 


1 ch 


B.2 Nomenclature 
Description 

Chamber calibration coefficient 

Nozzle diameter 

Wire diameter 

Chamber diameter 

Measured pressure differential 

Voltage differential across a Wheatstone bridge 

Gravitational constant 

Wire diameter Grashof number 

Average convective heat transfer coefficient 

Electrical current 
Fluid thermal conductivity 
Wire length 

Average Nusselt number 
Atmospheric pressure 
Static chamber pressure 

Prandtl number 

Heat transfer per unit length 

Electrical resistance 

Wire diameter Rayleigh number 

Chamber diameter Reynolds number 

Wire diameter Reynolds number 

4-Wire RTD reading 

Ambient temperature 

Chamber temperature 

Film temperature 

Wire temperature 


Units 

[] 

M 

M 

M 

in H 2 o] 

[V] 


[] 

Wlm 2 °K 

[A] 

[W/moK] 

M 

[] 

N / m 2 
N / m 2 

[] 

[W /m] 

[«] 

[] 

[] 

[] 

[*] 
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V 

Measured voltage 

[v] 

X 

Velocity magnitude (speed) 

[m / 5 ] 

Y 

Compressibility factor 

[] 

P 

Volumetric thermal expansion coefficient 

K~ l 

8 

Perturbation marker 

[] 

y 

Specific heat ratio 

[] 

A 

Variable resistance 

[n] 

P amb 

Ambient density 

kg / m 

P ch 

Chamber density 

kg / m 

P / 

Film temperature density 

kg / m 

Pw 

Wire temperature density 

kg / m 

CO 

Total probe resistance 

[«] 


B.3 Qualification of “Low-Speed” Anemometry 

“Typical hot wire probe calibration curves show a marked decrease in sensitivity at low velocities. 
This leads us to ask what the lower limit is for practical hot wire and hot film measurements” (Lomas, 5 
p. 66). We focus on answering this question in a concrete fashion. We know that if a probe is exposed to a 
quiescent fluid, natural convection currents develop and, as such, the probe loses thermal energy via 
convection in “still” air (from here out, we ignore conductive and radiant effects). When the probe is 
exposed to a low velocity flow, this natural convection cooling effect “masks” the forced convection 
cooling. We attempt to determine through analysis, and verify by calibration, when this masking becomes 
small enough that velocity can be properly determined. 

As a first cut at quantifying where a King’s law type relationship will break down, we cite a “rule of 
thumb” for measurements in air. In general, the minimum practical measurable velocity is in the 
neighborhood of 15 - 20 cm/s (Jorgenson 3 ). 

We next look at a relationship given by Collis and Williams 1 for large aspect ratio sensors ( L » D ). 


in in 1.85Gr w 


0.39 


\ 0.76 


T 

V 1 amb J 


The above relation describes a minimum Reynolds number for which a hot wire may be used. We define 
the Grashof number based on wire diameter as: 

Gr w = g$(T w - T amb )Dl /v 2 . 


The Grashof number is a “ratio of buoyancy to viscous forces” (Incropera and De Witt 2 p. 320). The 
familiar Reynolds number based on wire diameter is introduced as: 

Re w = xD w /v . 

Finally, we introduce (Incropera and DeWitt 2 p. 449): 

P _ _ \ P amb ~ Pw 
P w 1 amb ~ T\v 
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The above approximation is known as the volumetric thermal expansion coefficient. It is a 
thermodynamic property of the fluid that relates changes in density to changes in temperature at constant 
pressure. 

Suppose we are given the following representative data: 


T w = 525K 
T amb = 300 if 
T f = A\2.5K 

g = 9.81 m / s 
D = 5 pm 
L= 1mm 

v = 15.89x10 ~ 6 m 2 ls 
v f = 27.91xl0 _6 m 2 Is 
p w = 0.6647 kg! nr’ 

P amb= l - 16U kglm* 

= 0.8466 kg! nr’ 


The property values above are for air at STP (Incropera and DeWitt 2 p. 839). The probe dimensions are 
characteristic of the ones used in the lab. Combining the above definitions and solving for the minimum 
measurable fluid velocity, we find: 


x min 1-85 


A, 


U- 




V 


Pamb P ii 

P w T a mb — T\v J 


3 4 


( T w~ T amb)~o~ 


0.39 . 


,0.76 


T 

V 1 amb ) 


Numerically speaking, 


We have arrived at a contradiction to our “rule of thumb.” We now depart from accepting published 
suggestions and look at the problem from a more analytic point of view. Recall that the Reynolds number 
describes the relative importance of momentum effects to viscous effects. As mentioned earlier, the 
Grashof number is a ratio of buoyancy to viscous forces. This suggests that one may take the ratio: 

Gr w /Rel 

and attach some kind of significance to it. A moment of reflection suggests that the above ratio is a 
measure of the relative effect of natural convection (buoyancy) to forced convection. Thus, if the ratio is 
large, then natural convection effects are significant. Suppose we demand that: 

Gr D /Rej) « 1 . 

Given fluid properties and wire geometry, one can get a rough feel for a minimum measurable velocity. 
Turning back to the example at hand: 
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Gr w = 3.629 xl0“ 6 
Re;, = i 2 (0.0990 s 2 /m 2 ) 

If we boldly assume that natural convection is negligible provided that: 

Gr w /Re^ < 1/10 


we immediately solve for a lower limit of: 


i > 2 cm/s. 

Our analysis suggests that we may set another lower limit. 

For the sake of completeness, we include one more idea. We balance Nusselt numbers for forced and 
natural convection. Suppose the flow under consideration varies from 0 to 1 m/s. Evaluating the Reynolds 
number at the film temperature we find 


0<Re w , <0.18 

This range of Reynolds numbers is lower than what is commonly tabulated for a cylinder in cross 
flow. Nevertheless, using Reynolds analogy, one can extrapolate the behavior of Nusselt numbers from 
tabulated drag coefficient data. Looking at Incropera and DeWitt 2 (p. 343) one observes smooth, non- 
accelerating drag coefficient behavior for flow over a circular cylinder. With this in mind, we extrapolate 
on available Nusselt number data (Incropera and DeWitt 2 p. 344), and state: 

Nu w j orced = 0.989 Re®; 33 Pr 1/3 . 

For purely free convection we have: (Incropera and DeWitt 2 p.465): 

Nu w ^ee « CRa w n 


Where the Rayleigh number is defined as: 

Ra w = Gr w Pr . 

Evaluating the Rayleigh number in air at the film temperature we find: 

Ra w = 4.09x1 0 -7 . 

We have (Incropera and DeWitt 2 p. 465): 

Nu w free =0.675R« h °- 058 - 0.280. 

At this point we make a similar argument to what was done earlier. Suppose that we demand that: 

^ u w,free , 

« 1 . 

N u w, forced 

In other words we are demanding that the dominating mode of heat transfer is forced convection. Again 
we boldly suppose that forced convection dominates if: 
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Nu w jree = 0-675 Ra/ 058 
Nu w , forced 0.989 Re^,' 33 Pr 1/3 ~~ 

This assumption allows us to back out a minimum measurable velocity. Carrying out the algebra, we find: 

x > 1 1 cm Is. 

Given four different points of departure, we find four different lowest measurable velocities. 

However, all four are in the neighborhood of 0.1 mis. This gives one an idea of where to begin looking for 
non-linearities in a King’s law type calibration curve. In a different light, the last statement suggests that, 
given a data set, one can expect that all fluid velocities outside (and larger) than the neighborhood 0.1 mis 
can be made to lie on a straight line for the correct choice of the exponent in a King’s law curve. 


B.4 Low Speed Hot Wire Calibration 

Our in-house laboratory calibration facility consists of a micro-manometer, a RTD temperature 
sensing device, several voltmeters and a test chamber, as shown in figure B. 1 . The facility will accept any 
type of compressed gas. In the present case, we speak of bottled compressed air. Regarding figure B.2, the 
air enters on the left-hand side of the picture and passes over a baffle, which is immediately followed by 
an RTD (not shown). It then flows through a series of “quieting” screens. The air then passes through a 
large nozzle into a straight section (the “chamber”). It then passes through a second nozzle and exits into 
the room. 

The probe is inserted into the chamber and is aligned such that the wire is normal to the flow. Directly 
across from the probe is a static pressure tap, where the manometer is connected. A data point is collected 
by turning on the flow, measuring the resistance across the RTD, the voltage across the probe and the 
fluid level inside the manometer. Iterating this procedure for different flow rates produces a collection of 
data similar to those presented in table B. 1 . 

The next step is to convert the measured pressure to a velocity. According to Wilson, 9 the equations 
(which are specific to our calibration rig, but similar to commercial rigs) to make this conversion are: 
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r 

+ 3.7937 x 10 -15 

) 

The actual computation of velocity is tedious and will be omitted for the sake of brevity. However, 
the velocity-voltage relationship that results from transforming table B.l is interesting and is presented 
in figure B.3. We see a square root like behavior. If we choose the exponent, m, in King’s law to be 
m = 0.435, then we arrive at a result like figure B.4. Figure B.4 presents the result of several calibration 
runs. In this figure, one may observe that, upon casting the voltage and velocity data upon these 
coordinates, the results lie on a straight line. The linear approximation of the data can be expressed as: 

i 0435 =- 1 . 5931 + . 50 33 y 2 

Recall in our earlier analysis, we predicted that our minimum measurable velocity would be in the 
neighborhood of 0. 1 m/s. In other words, we would expect linear behavior into the neighborhood of 
0.1 m/s. The information presented in figure B.4 agrees with this prediction (for example, the left-hand 
end of the curve fit corresponds to a velocity of 4.8 cm/s). 

Regarding figure B.4, we note that for all abscissal values less than 4, the pressure differential is on 
the order of 0.0001 in. of water. This is problematic in the sense that the pressure measurement device is 
graduated by 0.001 in. of water. 

B.5 Calibration Error Analysis 

We have generated a calibration curve and, through it, we can relate voltage and velocity. We now 
ask; how good is our calibration relationship? This boils down to answering three questions; what is the 
error associated with the abscissa, the error associated with the ordinate and the error associated with the 
curve fit in figure B.4? 

We begin by determining the error associated with pressure-derived velocity measurements (the 
ordinate of fig. B.4). Following the methodology of Yavuzkurt, 10 we write (substituting our calibration 
device relationship for his): 


* clinch 
^ Ch y 


=1.1518-4.4131x10" 


-6 [ x ch^cli 


+ 1.7252x10" 


-1 1 f * clinch 


x ch 



We have incorporated the above equation into a Matlab code. 1 The user supplies the code with RTD 
and pressure measurements. The code solves the above polynomial equation (in x), strips off the 


'See appendix C. 
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unreasonable velocity values (the above equation is 3rd order, and, as such, has 3 zeros) and returns the 
chamber velocity. 

Several reduction equations are buried in the code. We state them now. 


T ch °C = 


RTDR- 200.43 
0.771429 


Where RTDR indicates the RTD reading observed on a digital voltmeter (where the DVM is assumed to 
be free of error). We state that a measured and converted resistance is within 0.05% of the actual 
temperature. We do not propagate any uncertainty in this calibration; instead, we make our range of 
potential room temperature fluctuations large enough to capture any RTD calibration errors. 

We also have: 


P ch - 


2 DelP 
4.019x10 


+ P 

_3 atm 


This is a conversion from the measured pressure to absolute chamber pressure. We assume that the 
conversion factor is free of error. 

Owing to the complexity of the above equation, we choose the following method for determining 
uncertainty associated with it (Strykowski 7 ): 

Step 1 

We calculate “average values” by feeding the above calibration equation the list of values shown in 
table B.2. We call the output of table B.2 the average velocity, x . 

Step 2 

We rewrite our calibration equation in the following form: 
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Observe that if we do not perturb any of the quantities (e.g., viscosity) in the above equation then we find 

x = x perturbed • 

Step 3 

We define the variable quantities: 

D ch =0.025654/?? ±0.000001m @ 95% confidence 
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d = 0.008128/77 ±0.000001/// @ 95% confidence 
3902/V/m 2 + 30( 

RTDR- 200 A3 


P atm = 9S902N /m 2 ± 3000N Im 2 @95% confidence 


T ch ~ 


Pch - 


.771429 
2(DelP±.0002) 


-±5°C @ 95% confidence 


+ P ntm @ 95% confidence 


3 1 


4.019x10” 

Implicitly linked to the chamber temperature is viscosity, which we state as: 

Vch( T ch) = 182.1xlO -7 IV o s/m 2 ±2.5x\0~ 7 N ° s/m 2 @ 95% confidence 


The chamber density can vary with both temperature and pressure through the relationship: 


r* / t* d \ _ Pch(Patin') 

Pch 1 1 ch >T C h > ~ “7 A, 

■*V7 ir-*- rh 


The only quantity that does not fall into the above framework is the chamber calibration coefficient: 


f xD ch 7 

+ 1.7252XKT 11 

( xD ch ^ 

l 

+ 3.7937 xlO -15 

( kD ch ' 

v ^ ch J 


V ^ ch J 


V ^ ch j 


C ch =1.1518-4.4131x10" 


The data provided by Wilson suggest that a safe uncertainty value for the chamber coefficient is: 

8 C ch = ±3% @ 95% confidence. 


We assume that the gas and gravitational constants are free of error. Lastly, a quick glance at a 
thermodynamics book will reveal that for all practical purposes, the specific heat ratio can be assumed to 
be constant over our range of temperature perturbation. 

Step 4 

We use the equation of step three with all things evaluated at their average value with the exception 
of the perturbed quantity. We calculate the deviation that results by perturbing each quantity 
independently (unless perturbing a quantity implicitly perturbs another i.e., pressure and density) and 
record their effects. 

Suppose, for the sake of clarity, we are perturbing temperature (but keep in mind that it could apply to 
any of the variable quantities). We make the following definitions: 


T ch + §T C h *T 

T, h - §T c h “ > X T ■ 


This is read as; perturbing temperature leads to a velocity perturbation based on temperature. The 
results that fall out of applying the equation of step 2 and the definitions of steps 3 and 4 are presented in 
table B.3. 

Step 5 

We define the total velocity perturbation due to temperature as: 
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1-4- 

+ 

1- 

x T 

X 



X 


2 


That is, we are taking the average of the magnitudes of the perturbations. Finally, we define the total 
effect of all perturbations as: 


Si 


pdv 



where i ranges over all of the variable quantities. The results are presented in table B.4. 

We have described the error associated with a pressure derived velocity measurement. We turn to 
uncertainties associated with fitting a line to a set of data points. This will add a second error, S Xf it . 

Turning back to our calibration data, we have the curve fit: 

i°- 435 = -1.5931 + 0.5033V 2 . 

According to Strykowski, 8 the error associated with a curve fit is described by the product of the 
standard estimate of error and a Student t tabular value at 95% confidence: 


SEEx x 


9,95% 


“ x 9,95% x | 


11 


Yj ( ' x i ~ x my 


(=i 


@95% confidence 


We transform the above linear relation back to x flt and compare this with our “average velocity,” x as 
presented in table B.2. This process gives the standard estimate of error as: 

SEE = 0.001437 


Looking to a Student t table, we find: 


x 9,95% = 2.262 

Thus, we expect the curve fit to differ from measured velocity be the amount: 
dx flt = SEE * x 9 95% = ±0.33% @ 95% confidence 

Flaving quantified the error associated with a pressure-derived velocity and the error associated with a 
curve fit, we now collect these two quantities and call it the total velocity error: 

/*2 ’2 

5 x tot (x) = -J5 xj it + 5x pdv @ 95% confidence 

The above result is presented in table B.5 and figure B.5. 

The last source of calibration error is due to voltage measurement uncertainty. Owing to the quality of 
our DVM, we take this error as negligible. 
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B.6 Further Errors 


We now mention some noteworthy sources of error that are beyond the scope of this appendix. 

Probe Misalignment 

It was assumed that the flow was normal to the wire. There are situations, both in calibration and in 
actual measurement where this may not be the case. Instead of measuring “w” only, now “v” and perhaps 
“w” are involved. There is no sure-fire way to deal with this type of error outside of using multi-wire 
probes in flows where more than one component of velocity is expected. 

Calibration Drift 

As time passes, the wire may collect some kind of filth. An accumulation of dirt influences the wire’s 
heat transfer coefficient and hence makes a calibration curve inaccurate. Moreover, owing to the fact that 
the wire is maintained at an elevated temperature, its thermophysical properties may change with time. In 
other words, its resistance-voltage-temperature relationship may change. How are these effects 
quantified? Who knows? The moral is that the wire should be calibrated often. 

Relative Humidity Effects 

The calibration analysis that was presented did not take this error into account. If significant 
variations are encountered, a correction should be made. For the conditions of the present experiment, this 
effect can be neglected. 

Changes in Ambient Pressure 

A change of ambient pressure between calibration time and measurement or a significant change over 
the course of an experiment can lead to errors. If these conditions are encountered, a correction should be 
made. For the period of this study, such corrections were considered to be unimportant. 

A final note: figure B.5 might lead the reader to believe that hot-wire velocity measurements are 
within 3% of the true value (except at low velocities). In reality, due to errors outside of calibration (e.g., 
the ones cited above and others), the error band is on the order of 5% (this value of 5% is typical of many 
references where hot-wire error in carefully conducted experiments is given). 
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B.7 Figures and Tables 



Figure B.l. — The hotwire calibration facility. 


5 groups of 80 mesh screens probe holders 
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0 435 Voltage (V) 



Velocity (m/s) 

Figure B.3. — Velocity vs. voltage as obtained from tables B.l and B.2. 



Figure B.4. — A typical calibration curve. 
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Figure B.5. — Percent error as a function of velocity at 95% confidence. 


TABLE B.I.— DIFFERENTIAL PRESSURE MEASUREMENTS 
(IN. OF WATER), VOLTAGE MEASUREMENTS (V) 
AND RTD RESISTANCE READINGS (OHMS) 


FPVO 

RTDR 

DelP 

1.9439 

217.36 

0.0004 

1.9892 

217.35 

0.0013 

2.018 

217.34 

0.0023 

2.0367 

217.33 

0.0035 

2.0887 

217.33 

0.0085 

2.1362 

217.34 

0.0176 

2.2067 

217.34 

0.0437 

2.2413 

217.34 

0.0651 

2.2757 

217.34 

0.0943 

2.3288 

217.3 

0.1609 

2.3832 

217.27 

0.2651 
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TABLE B.2.— CODE INPUTS AND OUTPUTS 


Code Input 


Code Output 

RTDR 

DelP 


xdot 

217.36 

0.0004 


0.067572 

217.35 

0.0013 


0.12177 

217.34 

0.0023 


0.16193 

217.33 

0.0035 


0.1997 

217.33 

0.0085 


0.31098 

217.34 

0.0176 


0.4471 

217.34 

0.0437 


0.70328 

217.34 

0.0651 


0.85745 

217.34 

0.0943 


1.0307 

217.3 

0.1609 


1.3432 

217.27 

0.2651 


1.719 


TABLE B.3.— PERTURBATION RESULTS 


D 


Atmospheric pressure 

Chamber pressure 

+ 

- 


+ 

- 

+ 

- 

0.067567 

0.067578 


0.066569 

0.068622 

0.082759 

0.047781 

0.12176 

0.12178 


0.11996 

0.12367 

0.13081 

0.11202 

0.16192 

0.16194 


0.15952 

0.16445 

0.16882 

0.15473 

0.19968 

0.19972 


0.19673 

0.20281 

0.20533 

0.19391 

0.31096 

0.31101 


0.30635 

0.31583 

0.31462 

0.3073 

0.44706 

0.44713 


0.44043 

0.45407 

0.44963 

0.44455 

0.70322 

0.70333 


0.69276 

0.71428 

0.70488 

0.70167 

0.85738 

0.85752 


0.84461 

0.87089 

0.85877 

0.85613 

1.0306 

1.0308 


1.0152 

1.0469 

1.0318 

1.0296 

1.3431 

1.3433 


1.3229 

1.3643 

1.344 

1.3423 

1.7189 

1.7192 


1.6931 

1.7462 

1.7197 

1.7184 





d 


Chamber constant 

Temperature and viscosity 

+ 

- 


+ 

- 

+ 

- 

0.067589 

0.067556 


0.0696 

0.065545 

0.068145 

0.066995 

0.1218 

0.12174 


0.12543 

0.11812 

0.12281 

0.12073 

0.16197 

0.16189 


0.16679 

0.15707 

0.1633 

0.16054 

0.19975 

0.19965 


0.20569 

0.19371 

0.2014 

0.19799 

0.31106 

0.31091 


0.32031 

0.30165 

0.31363 

0.30831 

0.44721 

0.44698 


0.46051 

0.43368 

0.45092 

0.44324 

0.70345 

0.7031 


0.72438 

0.68218 

0.70932 

0.69718 

0.85766 

0.85724 


0.88317 

0.83173 

0.86484 

0.84999 

1.031 

1.0305 


1.0616 

0.99979 

1.0396 

1.0217 

1.3435 

1.3428 


1.3835 

1.3029 

1.3549 

1.3314 

1.7195 

1.7186 


1.7706 

1.6675 

1.7341 

1.7038 
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TABLE B.4. — §i prfv =/(i) 


Velocity 

RSS % error 

0.067572 

26.1134943 

0.12177 

8.460908194 

0.16193 

5.5659778 

0.1997 

4.49661477 

0.31098 

3.666309389 

0.4471 

3.520029483 

0.70328 

3.484234107 

0.85745 

3.481403994 

1.0307 

3.481605379 

1.3432 

3.485216863 

1.719 

3.486751534 


TABLE B.5.— TOTAL ERROR ASSOCIATED WITH 
VELOCITY MEASUREMENT AT 95% CONFIDENCE 


Velocity 

RSS % error 

0.067572 

26.4434943 

0.12177 

8.790908194 

0.16193 

5.8959778 

0.1997 

4.82661477 

0.31098 

3.996309389 

0.4471 

3.850029483 

0.70328 

3.814234107 

0.85745 

3.811403994 

1.0307 

3.811605379 

1.3432 

3.815216863 

1.719 

3.816751534 
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Appendix C — Samples of C and Matlab Codes 

C.l Introduction 

This appendix serves to acquaint the reader with the coding routines used for data reduction. It will 
provide the curious reader with the necessary information to fully understand the data presented in this 
report. It will also be of use to the future researcher as a point of departure for new or modified coding 
routines. For reference, comments (lines ignored by the program) in C programs appear in the following 
form: /*comment*/. In Matlab programs comments take on the following form: %comment. 

C.2 A Simple Traverse Code in C 

/*This code allows the user to specify a position file that consists 
of absolute positions. In other words, you specify the locations that 
you want to the probe to go to. One thing to know is that the way the 
code is written, it takes the first point as a datum, i.e. the first 
point should be some initial point that can be identified (e.g. a 
wall, or a centerline or something) . Thus, before the code runs, move 
the probe to the initial point. As this code is written, one should be 
able to cut and paste in into some kind of shell/ scripting window and 
expect it to run.*/ 

/*The command to make the thing run is "gcc -o {executable filename) 
{codename} cib.o -lm. You must have the file "cib.o" in the same 
directory as the source code. */ 

/* The next two lines are library inclusions relating to the GPIB 
card. Depending on the machine that is being used, either may be 
appropriate. */ 

#include "/usr/local/include/sys/ugpib.h" 

/*#include<ugpib.h>*/ 

/* library inclusions*/ 

#include<math . h> 

#include<stdio.h> 

#include<unistd.h> 

#include<stdlib.h> 

#def ine BUFSIZE 65550 
#def ine STEPSIZE 400 
#def ine LOWSPEED 1000 

/* this is the number of points in the traverse, it must be changed 
as necessary*/ 

#def ine NUM_LOCATION 6 

/* beginning of code*/ 
main ( ) 

{ 


/* variable definitions */ 
int i; 
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int trav; 

int flag = 0; 

char command [ 20]; 

double xini = 0; 

double xcur [NUMJLOCATION] ; 

double xstep; 

double xabs; 

long step, wait; 

/* file pointer definitions */ 
FILE * fpposition; 

FILE * fpdataf ilename; 


/* Initialize Traverse for movement. You will have to change the 
traverse address as necessary. Also, you will have to change the 
traverse settings to suit your needs. Refer to the operators manual 
for this purpose. */ 

trav=ibf ind ( "dev3 " ) ; 
ibwrt (trav, "E , B" , 1 ) ; 
sprintf (command, " S1MLOWSPEED, ") ; 
ibwrt (trav, command, strlen (command) ) ; 

/* this tells the operator to move to the traverse to a datum. */ 
printf("move the traverse to zero displacement \n" ) ; 
printf("when this is the case, press l\n"); 

scanf ("%d", &flag) ; 

if ( f lag==l ) 

{ 

/* beginning of traverse movement code */ 

/* This directs the program to look for a position file. The 
position file should be in the same directory as the driving code 
(i.e. what you are looking at right now. */ 

/* The position file should have a name like "positionf ile" and 
should look like: /* 

/* 

0 

1 

2 

3 

3.1 

*/ 

/* you can specify any numbers you want. The first number need not 
be zero (we established a datum above) . You will need to count up the 
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number of positions in the file and insert this number in the 
definition "#define NUM_LOCATION {fill in the number here}" above, two 
notes: 1, all positioning is done in INCHES. 2. all positioning is 
ABSOLUTE (i.e. not incremental) */ 

fpposition=f open ("tester", "r") ; 

for (i = 0 ; i<NUM_LOCATION; i + +) 

{ 

f scanf (fpposition, "%lf ", &xcur [i] ) ; 

/*this reads in the absolute position file*/ 

} 


for (i = 0; i<NUM_LOCATION; i++) 

{0 

/* move to i'th point. */ 

xabs = xcur [ i ] ; 
xstep = xabs-xini; 

printf ( "moving from %lf to the new absolute position of %lf\n", 
xini , xabs ) ; 

printf ("step distance = %lf \n" , xstep); 


/* this just tests if the step size is valid */ 
if ( fabs (xstep) >= 0.0001) 

{ 

step = (long) (xstep*25 . 4*STEPSIZE) ; 

/* conversion to metric for the traverse */ 

if ( step < 0 ) 

/* moves the traverse in the negative direction */ 

{ 

sprintf (command, "C, IlM-%ld, R, C" , (step* (-1 ) ) ) ; 
wait=l-step/LOWSPEED; 

} 

if ( step > 0 ) 

/* moves traverse in positive direction */ 

{ 

sprintf (command, "C, IlM%ld, R, C" , step) ; 
wait=l+step/LOWSPEED; 

} 

ibwrt (trav, command, strlen (command) ) ; 
sleep (wait) ; 
sleep (10) ; 

/* you may have to play with the above sleep co"and. Essentially, it 
gives the traverse time to finish its move before the next command is 
issued. This isn't critical in this particular program, but in 
something else - like a motion/data ac program, it may be. For 
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example, you will want the traverse to be done moving before data 
collection begins */ 

} 

/* these next few commands are just so the user can observe the 
traverse movement at each point, if this code is incorporated into a 
larger one, more than likely, these lines will be removed */ 
printf("if everything looks OK press l\n"); 
scanf ("%d", Sflag) ; 

if (flag ! =1 ) 
exit ( 0 ) ; 

/*********************************************/ 
xini = xcur [ i ] ; 
f flush (stdout) ; 

/* end of traverse movement code */ 


C.3 A Simple Sampling Code in C 

/* This program collects hotwire data in raw form. */ 

/* Run and programmed by: DAA */ 

/*The command to make the thing run is "gcc -o {executable filename) 
{codename} cib.o -lm. You must have the file "cib.o" in the same 
directory as the source code. */ 

/* The next two lines are library inclusions relating to the GPIB 
card. Depending on the machine that is being used, either may be 
appropriate. */ 

#include<ugpib . h> 

/*#include" /usr/ local/ include/ sys/ugpib . h" * / 

/* standard inclusions */ 

#include<math . h> 

#include<stdio.h> 

#include<unistd.h> 

/* this defines "n" data points, you actually get n-1 useful points 
and a marker point at the end of the array. The marker is a "0" */ 

#def ine DATAPOINTS 65536 

/* beginning of Data Acq. */ 

main ( ) 

{ 

int j ; 
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int iotech, board, fluke ; 

/* gpib pointers to the iotech and interface board */ 

FILE * calfile; 

/* Where the calibration data is stored */ 
char * calfilename; 

/* The name of the calibration file */ 
static short data [DATAPOINTS] ; 

/* The waveform data, 256k data points */ 
static double voltage [ DATAPOINTS ] ; 

/* The voltage data, single channel */ 
double gain=1.0; 

printf ( "Acquiring data..."); 
calf ilename=" spectral . 4 " ; 

/* output filename */ 
calf ile=f open (calfilename, "a" ) ; 

/* Open data file for writing */ 

iotech=ibf ind ( "devl " ) ; 

/* find the iotech */ 
ibclr (iotech) ; 

/* Clear the iotech */ 
ibtmo (iotech, TIOOOs) ; 

/* Set timeout to 30 seconds */ 
ibwrt (iotech, "M4X" , 3) ; 

/* Clear buffer overrun mask if set */ 
ibwrt (iotech, "C1X" , 3) ; 

/* Use channel 1 on iotech */ 
ibwrt (iotech, "ROX" , 3) ; 

/* +/-1 volts range */ 
ibwrt (iotech, "G11X" , 4 ) ; 

/* Binary transfer mode */ 
ibwrt (iotech, "I5X" , 3) ; 

/* sample rate (2000 Hz) */ 
sleep ( 1 ) ; 

/* Give iotech a chance to catch up */ 

/* Acquire Data from Iotech, binary mode */ 
ibwrt (iotech, "NOX" , 3) ; 

/* Put iotech in FIFO mode */ 
ibwrt (iotech, "TOX" , 3) ; 

/* Trigger on talk */ 
sleep ( 1 ) ; 

printf ( "Acquiring data..."); 
f flush (stdout) ; 

ibrd (iotech, (char *)data, 2*DATAPOINTS) ; 

/* 2 bytes per data point, 1 channels, 256k data points */ 

printf ( "done . (%d bytes collected) \007\n" , ibcnt) ; 

f flush (stdout) ; 
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for ( j =1 ; j <=DATAPOINTS ; j ++ ) 

{ 


voltage [ j ] = ( (2.0/ gain) /60000.0) *data [ j ] ; 

/*the "2" is the voltage range, the "60000.0" is a conversion from 
binary*/ 

fprintf (calfile,"%lf \t\n", voltage [ j ] ) ; 

/*write voltage into a datafile*/ 

} 

fclose (calfile) ; 

} 


C.4 Low Speed Hot-Wire Calibration Code in Matlab 

%This pertains to the hotwire calibration facility in lab 281. The 
%user of this code must carefully follow these instructions to obtain 
%proper results. This code pertains to the "chamber" of the 
%calibration device as defined in Wilson's thesis (see Appendix B) . 

%This script file will take excel hotwire calibration data from a 
user %specified data file and perform the operations outlined is 
Wilson's %thesis to convert measured temperature and pressure readings 
to %velocity data, the only reason that Matlab has been used is 
because %the governing equations for flow within the calibration 
device 

%require finding roots of a 3rd order polynomial, 2 of the roots are 
%incorrect the third root is the desired velocity. The difference 
%between the desired velocity and the two false ones is obvious upon 

%inspection . 


TP = xlsread ( ' Book3 ' ) ; 

%above command gets the data from excel. Make sure to format excel 
%such that temperature sensor readings are in column 1 of the passed 
%array and pressure readings are in column 2. (temp to the left of 
%press.) Insert your excel filename between the single quotes. 

len = length (TP) ; 

%Determines the number of entries from the excel column vector. 

T = TP (1 : len, 1) ; 

%takes all the temperature readings and stores them in a column 
vector 

P = TP (1 : len, 2) ; 

%takes all the pressure readings from TP and puts them in a pressure 
column vector 

clear TP 

% clears out TP 
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%user defined constants 
P_atm = 97764.0; 

%the user should check this at calibration time 
gamma = 1.4; 

%specific heat ratio 
D_CH = .0256; 

%chamber diameter 
d_noz = .008128; 

%nozzle diameter 
R_air = 283; 

%metric gas constant for air 
RTD_coeff = .799178; 

%slope of the RTD calibration curve 
RTD_baseline = 200.49; 

%intercept of the RTD calibration curve 

%you must determine the voltage/temp relation between whatever type 
of sensor you are using 
mu = 184 . 6*10 A -7; 

%viscosity of air 


P_CH = 2 * P / (4 . 015*10 A -3) + P_atm; 

%A column vector. The above line converts the directly read pressure 
%measurement (in inches of water) to Pascals and adds on atmospheric. 

%Notice the factor of 2, this is because a U-tube manometer was 
used . 

T_CH = (T - RTD_baseline) /RTD_coef f + 273.15; 

%A column vector. This converts from RTD readings to temperature in 
%Kelvin 

for i = 1 : len 

%if you get a divide by zero message on the front end, you may have 
to %change the index 


rho_CH = P_CH (i) / (R_air*T_CH (i) ) ; 
%point by point determined constant 

sigma = (d_noz/D_CH) A 2 ; 

%constant 

omega = (P_atm/P_CH (i) ) A (1/gamma) ; 
%point by point constant 
chi = rho_CH*D_CH/mu ; 

%point by point constant 


alpha = 1.1518 * sigma * omega; 

beta = -4 . 4131*10 A -6 * chi * sigma * omega; 


NAS A/CR— 2006-2 14131 


177 



delta = 1 . 7252*10 A -11 * chi A 2 * sigma * omega; 
phi = 3 . 7937*10 A -15 * chi A 3 * sigma * omega; 

%point by point constants 

epsilon = ( ... 

(2 * (P CH(i)-P_atm) * gamma * ( 1 - (P_atm/P_CH (i) ) A ( (gamma - 

1 ) /gamma) ) ) . . . 

/. . . 

( rho_CH * ( 1 - ( P_atm / P_CH(i) ) A (2/gamma) * (d_noz/D_CH) A 4 ) 

* (gamma - 1) * (1 - (P_atm/P_CH (i) ) ) )... 

) • • • 

A (1/2) ; 

%point by point determined constant 


polynomial = [phi*epsilon delta*epsilon (beta*epsilon-l ) 
alpha*epsilon] ; 

polyroots ( i , 1 : 3 ) = transpose (roots (polynomial )) ; 


end 


dlmwrite ('Data.xls', poly roots, ' \ t ' ) 

%this is the output to an excel file named "Data". There will be 3 
%columns in the output. Two of them will be obviously incorrect values 
%of velocity, the third one is the one to keep. 

C.5 The Oscillatory Hot-Wire Traverse Code in C 

/* This program is written to take velocity measurements in the disc 
and tube spaces. Date: 04-12-2003. As arguments, it takes a list of 
absolute positions. These should be presented as a column vector of 
positions. The program will move the probe to the position values that 
you specify. The program also asks for the slope and intercept of a 
kings law calibration curve and the ambient temperature at calibration 
time. These values are input at the command line */ 

/* to run this code type gcc -o {executable filename} {code 
filename} cib.o -lm at the command line. The file cib.o must be in the 
same directory as the code */ 

#include "/usr/local/include/sys/ugpib.h" 

/*#include<ugpib.h>*/ 

/* select either of the above definitions depending on which machine 
you are using */ 

#include<math . h> 

#include<stdio.h> 

#include<unistd.h> 

#include<stdlib.h> 
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#def ine BUFSIZE 65550 
#def ine STEPSIZE 400 
#def ine LOWSPEED 1000 

/* CHANGE AS NECESSARY */ 

#def ine NUM_LOCATION 73 

/* count up the number of positions in the position file and stick 
it in here */ 

/* CHANGE AS NECESSARY to account for rotational frequency*/ 

#def ine NUM_CYCLE 102 

/* the number of cycles to be averaged */ 

#def ine NUM_POINT 1000 

/*30 RPM -- change IOTech to I7X*/ 

/*#def ine NUM_POINT 857*/ 

/* 7 0 RPM -- change IOTech to I6X*/ 

/* The number of points per one cycle, must be changed to 
accommodate for rotational frequency. We adjust NUM_POINT/sampling 
frequency such that it is identically equal to 1/rotational frequency. 
For example, 30 rpm => . 5 Hz rotational frequency => 1/rf = 2 === 
NUM_POINTS/sampling frequency */ 

main ( ) 

{ 

int i , j , k , 1 ; 
int flag = 0; 
int iotechl , iotech2 ; 
int trav, sequence ; 
char command [20]; 

FILE * fpposition; 

FILE * fpvel; 

FILE * fpdataf ilename; 

/* traverse variables */ 
double xcur [NUM_LOCATION] ; 
double xini = 0; 
long step, wait; 
double xstep; 
double xabs; 

double m,b; 

/*m = slope, b = intercept */ 
static double tdc_voltage; 

static short voltl [NUM_POINT] , volt2 [ 2*NUM_POINT] ; 
static double tdc [NUM_POINT] , pulse [NUM_POINT] ; 
static double vel_volt [NUM_POINT] ; 
static double vel [NUM_POINT] ; 

static double sum_vel [NUM_POINT ] , ave_vel [NUM_LOCATION] [NUM_POINT] ; 
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double gain=1.0; 

/*must change as is necessary*/ 
double Tsensor=250; 

/*this may not be right depending on overheat ratio */ 
double Tfluid; 

/* amb . air temp measured before the run */ 
double Tcalib; 

/*temperature at calibration time*/ 

/* Create filename*/ 
char filename [50 ] ; 

/* Get the calibration constant for the hot-wire*/ 
printf ( "enter the slope associated with the king's law calibration 
curve . \n" ) ; 

scanf ("%lf", &m) ; 

printf ( "enter the intercept associated with the king's law 
calibration curve\n") ; 
scanf ( "%lf " , &b) ; 

printf ( "enter the ambient temperature at calibration time\n"); 
scanf ( "%lf " , &Tcalib) ; 

printf ( "enter the current room ambient temperature\n" ) ; 
scanf ("%lf", &Tf luid) ; 

/* Initialize IOTECH1 for velocity */ 
iotechl=ibf ind ( "devl" ) ; 

/* change device number accordingly */ 
ibclr ( iotechl ) ; 

/* Clear the iotech */ 
ibtmo ( iotechl , T10 Os ) ; 

/* Set timeout to 100 seconds */ 
ibwrt ( iotechl , "M4X" , 3 ) ; 

/* Clear buffer overrun mask if set */ 
ibwrt ( iotechl , "C1X" , 3 ) ; 

/* Use channel 1 */ 
ibwrt ( iotechl , "R2X" , 3 ) ; 

/* +/-5 volts range */ 
ibwrt (iotechl, "G11X", 4) ; 

/* Binary transfer mode, little-endian */ 
ibwrt ( iotechl ," I7X" , 3 ) ; 

/* sample rate 7 => 500 Hz */ 
sleep ( 1 ) ; 

/* Give iotech a chance to catch up */ 

/*Initialize IOTECH2 for TDC and pulse*/ 
iotech2=ibf ind ( "devl4 " ) ; 
ibclr ( iotech2 ) ; 

/* Clear the iotech */ 
ibtmo ( iotech2 , T10 0s ) ; 

/* Set timeout to 100 seconds */ 
ibwrt ( iotech2 , "M4X" , 3 ) ; 
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/* Clear buffer overrun mask if set */ 
ibwrt (iotech2 , "C1X" , 3 ) ; 

/* Use channels 1 */ 
ibwrt (iotech2 , "R2X" , 3 ) ; 

/* +/-5 volts range */ 
ibwrt (iotech2,"GHX",4) ; 

/* Binary transfer mode, little-endian */ 
ibwrt (iotech2 , " I7X" , 3 ) ; 

/* sample rate */ 
sleep ( 1 ) ; 


printf("move the traverse to zero displacement^" ) ; 
printf("when this is the case, press l\n"); 
scanf ("%d", Sflag) ; 
if (flag ! =1 ) 
exit ( 0 ) ; 


/* Initialize Traverse */ 
trav=ibf ind ( "dev3 " ) ; 

/* change as is necessary */ 
ibwrt ( trav, "E , B" , 1 ) ; 
sprintf (command, " S1ML0WSPEED, ") ; 
ibwrt (trav, command, strlen (command) ) ; 


/* CHANGE AS NECESSARY - this is the position file */ 
fpposition=f open ( "tubelocationf ile" , "r" ) ; 

/* CHANGE AS NECESSARY - this is the position file */ 


for (i=0 ; i<NUM_LOCATION; i++) 

{ 

f scanf (fpposition, "%lf ", &xcur [i] ) ; 

/*this reads in the absolute position file update NUM LOCATION to 
match the file*/ 

printf ( "%lf \n" , xcur [i ] ) ; 

} 


for (i=0; i<NUM_LOCATION; i++) 

{ 

/*move to ith point.*/ 

xabs = xcur [i ] ; 
xstep = xabs-xini; 

printf ( "moving from %lf to the new absolute position of %lf\n", 
xini , xabs ) ; 

printf ("step distance = %lf \n" , xstep); 
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if ( fabs (xstep) >= 0.0001) 

/* this just tests if the step size is valid */ 

{ 

step = (long) (xstep*25 . 4*STEPSIZE) ; 

/* conversion to metric for the traverse */ 

if ( step < 0 ) 

/* moves the traverse in the negative direction */ 

{ 

sprintf (command, "C, IlM-%ld, R, C" , (step* (-1 ) ) ) ; 
wait=l-step/LOWSPEED; 

} 

if ( step > 0 ) 

/* moves traverse in positive direction */ 

{ 

sprintf (command, "C, IlM%ld, R, C" , step) ; 
wait=l+step/LOWSPEED; 

} 

ibwrt (trav, command, strlen (command) ) ; 
sleep (wait) ; 
sleep (10) ; 

} 

xini = xcur [ i ] ; 
f flush (stdout) ; 


/* Open a datafile for writing */ 

fpdataf ilename=f open ( "dataf ilename" , "w") ; 


/* CHANGE AS NECESSARY - This is where the data is written to */ 

fprintf (fpdataf ilename, "/h/adolf son/research/osc_work/2 . 125- 
30/tube/%lf\n" , xcur [i] ) ; 

/* CHANGE AS NECESSARY - This is where the data is written to */ 
fclose (fpdatafilename) ; 

fpdataf ilename=f open ( "dataf ilename" , "r") ; 
fscanf (fpdatafilename, "%s", filename) ; 
printf("Data file name is %s \n" , filename) ; 
fclose (fpdatafilename) ; 

fpvel=fopen (filename, "a") ; 

for (1=0; 1<NUM POINT; 1++) 
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sum vel [ 1 ] =0 . 0 ; 


/* Trigger IOTECH data acquisition*/ 
ibwrt (iotechl,"N0X",3) ; 

/* Put iotech in FIFO mode */ 
ibwrt (iotechl,"T0X",3) ; 

/* Trigger on TALK */ 

ibwrt (iotech2,"N0X",3) ; 

/* Put iotech in FIFO mode */ 
ibwrt (iotech2,"T0X",3) ; 

/* Trigger on TALK */ 
sleep ( 1 ) ; 

for ( j =0 ; j <NUM_CYCLE ; j ++ ) 

{ 

printf (" the %d th cycle \n",j); 


print f ( "Acquiring data. . . \n" ) ; f flush ( stdout ) ; 
ibrd ( iotechl , (char *) voltl, 2*NUM_POINT) ; 
printf ( "done . (%d bytes collected by IOTECH 

1) \007\n", ibcnt) ; f flush ( stdout ) ; 

ibrd ( iotech2 , (char *) volt2, 2*NUM_POINT) ; 
printf ( "done . (%d bytes collected by IOTECH 

2) \007\n", ibcnt) ;fflush( stdout ) ; 

for (k=0 ; k<NUM_POINT; k++) 

{ 

vel_volt [ k] = ( (10 . 0/gain) / 60 000 . 0) * voltl [k] ; 
tdc [k] =volt2 [ k] *(10. 0/60000); 

} 

for (k=0; k<NUM_POINT; k++) 

{ 

/* this converts voltage to velocity */ 

vel [ k] =pow ( (b+m*vel_volt [k] *vel_volt [k] ) ,2.29885) ; 
vel [k] =vel [k] *pow( (Tsensor-Tcalib) / (Tsensor-Tfluid) , 0.5) 
f printf ( fpvel ,"%lf\t%lf\n " , tdc [ k] , vel [k] ) ; 


printf ( " \n" ) ; 

} 

f printf ( fpvel , " \n" ) ; 
f close ( fpvel ) ; 

} 

fclose (fpposition) ; 

} 
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C.6 A Reduction Code for Oscillatory Flow in C 

/* this code takes the data from The Oscillatory Hot-Wire Traverse 
Code and ensemble averages it. Inputs to this code are the hot-wire 
derived velocity data files from the aforementioned code and a list of 
positions. It is best to number the velocity files to correspond with 
the position files. The reason being is that when the code accepts the 
list of positions, it reads and stores it from top to bottom. Thus, 
when the code reads "position 1 " from the position file, it goes 
looking for a velocity output file with the name "position 1", when it 
reads "position 2 " it goes looking for a file of the same name, etc. 

*/ 


/* To run this code, type gcc -o {executable filename} {code 
filename} -lm at the command line */ 

/* standard inclusions */ 

#include<math . h> 

#include<stdio.h> 

#include<stdlib.h> 

#def ine NUM_X 102000 

/* defined as 102 cycles* ( sampling f requency*rotational period)*/ 
#def ine NUM_CYCLE 80 

/* the number of cycles to be ensemble averaged. The code will read 
from top to bottom until the number of cycles has been reached */ 

#def ine NUM_FILE 73 

/* number of sampling locations in traverse */ 


main ( ) 

{ 

int i, j , k, N, M, 1 ; 

/*counter variables*/ 
int initial , index; 
int Npoints; 

FILE * fpold; 

/* this is the raw hotwire-TDC sensor data */ 

FILE * fpnewV; 

/* this is the velocity-crank angle output */ 

FILE * fpRMS; 

/* this is the Tl-crank angle output */ 

FILE * fpVf ilename; 

/* this will eventually be a pointer to a list of filenames in a 
traverse */ 


double data; 
double xl [NUM_X] ; 

/* an array of signals from the TDC sensor */ 
double tdc[NUM X] ; 
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double pulse [NUM_X] ; 
double Vel [NUM_X] ; 

/* an array of velocity signals */ 
double f[NUM_X]; 
double sumV,sumRMS; 

/*intermediate sums for average/TI calcualtions */ 
double aveV [ 720 ] , aveRMS [ 720 ] ; 

/* arrays of average values and TI ' s */ 
double location [NUM_FILE] ; 

char dataf ilename [50] ; 


/*******************************************************/ 

fpnewV=fopen ( "80cycVelDataCooler8-9" , "a" ) ; 

/* this is the output velocity file */ 
fpRMS=f open ( " 80cycVelRMSCooler8-9 " , "a" ) ; 

/* this is the output TI file */ 

for(i=0;i<720;i++) 

{ 

fprintf (fpnewV, "%d\t" , i) ; 

/* this makes 720 points in the output files i.e. half degrees */ 
fprintf ( fpRMS , "%d\t " , i ) ; 

} 

fprintf (fpnewV, "%d\n" , i) ; 

/* puts a newline between the angle index and the velocity row */ 
fprintf ( fpRMS , "%d\n" , i ) ; 

/* This is where the list of positions is read in */ 

fpVf ilename=f open ( " listof positions" , "r" ) ; 
for (i=0; i<NUM_FILE; i++) 

{ 

fscanf ( fpVf ilename, "%lf\n", &location [i] ) ; 
printf ( "%lf \n" , location [ i ]) ; 

} 


for (1=0; 1<NUM_FILE; 1++) 

{ 


sprintf (dataf ilename, "%lf",location[l] ) ; 
fpold=fopen (dataf ilename, "r") ; 

/*fpold=fopen ("VelCO . 0 0 0 0 0 0 " , "r " ) ; */ 

printf ("The Data being processed was sampled at a distance of %lf 
inches from the target impingement plate \n" , location [1 ]) ; 
printf ("Processing. . . ., please wait ! \n" ) ; 
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for ( i=0 ; i<NUM X;i++) 


{ 

f scanf ( fpold, "%lf\t%lf\n", &xl [i] , &Vel [ i ] ) ; 

} 


/* printf ("%lf\t%lf\t%lf\t%lf\n", xl [ 0 ] , x2 [ 0 ] , x3 [ 0 ] , x4 [ 0 ] ) ; */ 
/*f scanf (fpold, "%lf\t%lf\t%lf\%lf\n", &xl [i] , &x2, &x3, &x4) ;*/ 

i = 0; 

while (xl [ i ] >0 . 5 ) 

{ 

i=i + l ; 

/*f scanf (fpold, "%lf\t%lf\t%lf\%lf\n", &xl [i] , &x2, &x3, &x4) ; */ 

} 


N=i ; 

/*printf ("%d\t%lf\t%lf\n",N,xl [N] ,x2 [N] ) ;*/ 
index=0 ; 

for ( j =0 ; j <NUM_CYCLE ; j ++ ) 

{ 

Npoints=0 ; 
initial=index; 


while((xl[index+N]<l&&xl[N+index+l]>l) | | (xl[index+N]>l&&xl[N+index+l]> 
1 ) | | (xl [ index+N] <1 &&xl [N + index+1 ] <1 ) ) 

{ 

tdc [index] = j ; 

/*variableR [index] =x2 [index+N] ; 
variableC [ index] =x3 [ index+N] ; 
variableH [ index] =x4 [index+N] ; * / 


/*printf ("%lf\t%lf\t%lf\n", tdc [index] , TempR [index] , TempC [index] ) ; */ 
index=index+l ; 

Npoints=Npoints+l ; 


tdc [index] =j ; 

/*variableR [index] =x2 [index+N] ; 
variableC [ index] =x3 [ index+N] ; 
variableH [ index] =x 4 [ index+N] ; * / 
Npoints=Npoints+l ; 
index=index+l ; 

/ *printf (" index is %d\n" , index); 
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printf ( "Npoints is %d\n" , Npoints ) ; * / 


for (k=0 ; k<Npoints ; k++ ) 

{ 

if(((k+l)*720. O/Npoints-O .5) >floor ( (k+1) *720 . 0 /Npoints) ) 

{ 

pulse [k+initial]=floor( (k+1) *720. 0 /Npoints ) + 1 ; 

} 

else 

pulse[k+initial]=floor( (k+1) *720.0 /Npoints ) ; 


f printf (fpnewV, "%lf\t", location) ; 

/* calculate ensemble averaged velocities */ 
for(i=0;i<720;i++) 

{ 

M=0 ; 

sumV=0 . 0 ; 

for (k=0 ; kcindex; k++) 

{ 

if (fabs (pulse [k] - (i+1) ) <0 . 1) 

{ 

sumV=sumV+Vel [k+N] ; 

M=M+1; 

} 

} 

aveV [ i ] =sumV/M; 

f printf ( fpnewV, "%lf \t " , aveV [ i ] ) ; 

} 

/*Calculate TI levels*/ 

f printf (fpRMS, "%lf\t", location) ; 

for (k=0 ; k<index; k++) 

{ 

for (i = 0;i<720;i++) 

{ 

if (fabs (pulse [ k] - (i+1) ) <0 . 1 ) 

{ 

f [ k] =Vel [ k+N] -aveV [ i ] ; 

} 



for (i = 0;i<720;i++) 

{ 

M=0 ; 
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sumRMS=0 . 0; 

for (k=0 ; k<index; k++) 


if (fabs (pulse [k] - (i + 1) ) <0 . 1) 

{ 

sumRMS=sumRMS+pow ( f [ k] , 2 ) ; 
M=M+ 1 ; 


aveRMS [i] =sqrt (sumRMS/ (M-l) ) /a veV [i] ; 
/*this isn't RMS, it is TI*/ 

fprintf ( fpRMS , " %lf \t" , aveRMS [ i ] ) ; 

} 


fprintf ( fpnewV, " \n" ) ; 
fprintf ( fpRMS , " \n" ) ; 
f close ( fpold) ; 


} 

f close (fpnewV) ; 
f close ( fpRMS ) ; 

/*fclose ( fpVf ilename) ; */ 

} 


C.7 A Fourier Transform Code with Power-Spectrum Plotting in Matlab 

%This code takes an array of sampled data from excel, performs a 
discreet Fourier transform on it and packages the result in terms of 
power and frequency. Also provided is a Matlab power/spectral plot. 

filename = input ('enter your excel filename exactly as it appears 
(without the .xls) : ','s'); 

signal = xlsread (filename) ; 

%The above commands get the data from excel. Make sure to format 
excel %such that input data is in column vector form. Insert your 
excel %filename between the single quotes 

N = length ( signal ) ; 

%number of sample points for a N point DFT 

fs = input ('enter the data sampling frequency (in Hz).\n'); 

%this asks for and stores the sampling frequency. 

b = 1 : N; 

% bins, necessary for DFT computation and power spectra plotting 
% looks like: b = [1 2 3 . . . N] 

Ts = 1/fs; 

% sampling interval in seconds 
ts = Ts* (b-1 ) ; 

% discreet time points 
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% looks like: ts = Ts* [0 12 3 ...N-l] 

%computation of FFT and power spectra 
SIGNAL = fft (signal) ; 

%DFT of x 

pwr = SIGNAL. *conj (SIGNAL) /N; 

%this computes the "power" of the signal 
f rs = (b-1 ) /N* f s ; 

%plot of the results 

%this presents a space/time plot of the input signal. 

subplot (2,1,1) ; 

plot (ts, signal) 

xlabel ( ' Time , s') 

ylabel ( ' time domain signal ' ) 

grid on 


%this presents a power/frequency plot 
subplot (2,1,2) ; 
loglogffrs, pwr) ; 
grid on; 

xlabel ('Frequency, Hz'); 
ylabel ( ' Power ' ) ; 

%this is a verification of Parseval ' s theorem. The theorem states 
that %the sum square of the Fourier coefficients should equal the 
average %value of the signal squared over the time length of the 
sample. Here %is that check, for further details see any text on 
Fourier analysis. 

disp('this is the verification of Parsevals theorem'); 
format long 

[sum(pwr) norm ( signal ) A 2 ] 
format 

disp('the above numbers should be identical up to machine accuracy, 
if they do not then something is wrong and the results should be 
discarded ' ) ; 


C.8 A 3-D Plotting Code in Matlab 

%%%%%%%%%%this gets a matrix of values of crank angle and average 
%%%%%%%%%%velocities or TI's. It looks like: 

o. 2- 2- 2- 2- 2- 2- 2- 2- 2- D 1 o Q A i o n 

ooooooooooU -L Z. O e -± . . . / Z U 

S- S- S- S- S- S- S- S- S- S-tt-I t tO .7-790 

ooooooooooV_L VZ ... V/ZU 

2'2'2'2'2'2'2'2'2'S' 

oooooooooo ... 


%in short, this code will turn the results of the code "A Reduction 
Code for Oscillatory Flow in C" and make a pretty picture out of it. 
This particular code is specific to a tube profile, but following the 
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comments, a user will have little difficulty switching to a disc 
profile . 

%geometrical constant - this is specific to the tube 
D_tube = 8.465; 

%this scans in the position file and figures out how many spatial 
%positions there were 

filename = input ('enter position filename exactly as it appears : 

' , ' s ' ) ; 

posfile = dlmread (filename, ' \t ') ; 
clear filename; 
numpos=length (posfile) ; 

%here we take the tube or disc profile and apply the proper 
normalizing %processes and, in the case of a disc profile, a position 
shift (due to %the thickness of the probe) . We also scan in the file 
(or files in %the disc case) where data originates from. 

posfile = posf ile+D_tube/2 ; 

%this makes the centerline the zero position (dia->rad) 
posfile = posf ile/D_tube; 

%normalization r->r/D 

%this scans in the velocity - crank angle output file 

filename = input ('enter velocity output filename exactly as it 
appears : ' , ' s ' ) ; 

signal = dlmread (filename, ' \t ') ; 
clear filename; 


%this scans in the TI - crank angle output file 

filename = input ('enter TI output filename exactly as it appears 

' , ' s ' ) ; 

TI_signal = dlmread ( filename ,' \t ') ; 
clear filename; 

%manipulations on the velocity and TI matrices. Essentially, we 
%remove the column corresponding to count = crank angle = 0 because it 
%is garbage, we also strip off the l->720 numbering sequence and keep 
%in mind that as you step rightward in the matrix, you are increasing 
%in index (or crank angle) . 

vel_ca = signal (2 : numpos+1 , 2 : 721 ) ; 

clear signal; 

maxvect = max(vel_ca) ; 

max_vel = max (maxvect) ; 

clear maxvect 

TI_ca = TI_signal (2 : numpos+1 , 2 : 721 ) ; 
clear TI_signal; 
maxvect = max(TI_ca); 
max TI = max (maxvect); 
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clear maxvect 


%this chops off the ends of the data set to eliminate wall 
conduction %effects and blowing through the access hole 
x = 12; 

chop_vel = vel ca (2 : numpos-x, 1 : 720 ) ; 
position = posfile (2 : numpos-x) ; 
crankangle = 0.5:0.5:360; 
maxvect = max ( chop_vel ) ; 
max_vel = max (maxvect) ; 
clear maxvect 

%this makes a surface plot of velocity vs. radial position and crank 
%angle 

crankangle = 0.5:0.5:360; 

surf (crankangle, position, chop_vel) 

axis([0.5,360,min(posfile) , max (posfile) , 0 , max_vel ] ) 
ylabel ( ' r/D ' , 'FontSize',20) 
xlabel ( ' \ theta ' , ' Fonts ize ' , 20) 
z label (' | | u | | ', 'FontSize',20) 


%this makes a TI contour map 

%[X,Y] = meshgrid (crankangle, position) ; 

%[C,h] = contour (X, Y, chop_TI ) ; 

%h = clabel (C, h, ' FontSize ' , 14) ; 

%set (h, ' BackgroundColor ' , [ 1 1 . 6 ] , ' Edgecolor ' , [ . 7 .7 .7]) 

%xlabel ( ' \ theta ', 'FontSize',20) 

% ylabel ( ' r/D ' , ' Fonts ize ' , 20) 

C.9 A Code to Determine Wall Shear Stress as a Function 
of Crank Angle in an Oscillatory Flow in Matlab 

%given an array of velocity values like as delivered from the code 
of %A3 . 6 : 

%0 1 2 3 4 ... 720 
%vl v2 ... v720 

Q_ 

O ... 

%this program solves for wall shear stress and the error associated 
%with it as well as plots the results 

%geometrical constant - change as necessary 
disc_spacing = 0.127; %5 inch 

%this scans in the position file and figures out how many spatial 
positions there were it also takes the positions (in inches) to 
meters . 

filename = input ('enter position filename exactly as it appears : 

' , ' s ' ) ; 

posfile = dlmread (filename, ' \t ') ; 
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posfile = 0 . 0254*posf ile; 
clear filename; 
numpos=length (posfile) ; 


%this scans in the velocity - crank angle output file 
filename = input ('enter velocity output filename exactly as it 
appears : ' , ' s ' ) ; 

signal = dlmread (filename, ' \t ') ; 
clear filename; 

%manipulations on the velocity matrix. Essentially, we remove the 
%column corresponding to count = crank angle = 0 because it is 
garbage . 

%we also strip off the l->720 numbering sequence and keep in mind 
that %as you step rightward in the matrix, you are increasing in index 
(or %crank angle) . 

vel_ca = signal (2 : numpos+1 , 2 : 721 ) ; 
clear signal; 

%here we find the number of the index associated with the axial 
%location x/S = 0.100 

i=l; 

while (posfile (i) /disc_spacing) <0.1 
i=i+l ; 
end 

%now we shorten up the position file 
posfile = posf ile (1 : i) ; 

%now we eliminate far wall velocity data 
vel_ca = vel_ca (1 : i, 1 : 720) ; 

for j=l : 720 

%the following line finds the minimum velocity at each crank angle 
[C ( j ) , I ( j ) ] =min ( vel_ca (1 : i, j ) ) ; 

%now we take the spatial location of that minimum value and the next 
5 %locations 

y ( 1 : 5 , j ) = posfile (I (j ) : I ( j) +4) ; 

%now we take the corresponding velocities 
u ( 1 : 5 , j ) = vel_ca (I ( j ) : I ( j ) +4, j ) ; 
end 

%Average velocity and positions 
yavg = sum(y)/5; 
uavg = sum(u)/5; 

%place holders 
Y (1 : 720, 1 : 2) = 0; 

U (1 : 720, 1 : 2) = 0; 
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%fixed origin and average velocity/position in two arrays to 
facilitate %next step 
for i = 1:720 

Y ( i , 2 ) = yavg (i); 

U ( i , 2 ) = uavg ( i ) ; 
end 

%linear fit of average velocity to average position 
for i =1:720 

p (i, 1 : 2) = polyf it (Y (i, 1 : 2 ) , U (i, 1 : 2) , 1) ; 
end 

mu = 184.6*10 A -7; % units: N*s/m A 2 

%computation of wall shear stress 
WSS = mu*p (1 : 720, 1) ; 
theta = . 5 : . 5 : 360 ; 

%now for the error check. Now that we have wall shear stress as a 
%function of crank angle, we can recast our velocity on U+, Y+ 
%coordinates . That is, coordinates that are specific to the near-wall. 
%Then, we take the converted velocity data and compare it to a line 
(in %near wall coordinates) with a slope of 1. This should give us an 
idea %of the error associated with the original line-slope 
determination of %shear stress made above. 

%density 
rho = 1.1614; 

%shear velocity 

for i = 1:720 

vs tar ( i ) = (WSS ( i ) /rho ) A . 5 ; 

end 

%viscosity 
nu = mu/ rho; 

%conversion to u+ y+ co-ordinates 
for k = 1 : 720 

yplus(l:5,k) = y ( 1 : 5 , k) *vstar ( k) /nu ; 
uplus(l:5,k) = u ( 1 : 5 , k) /vstar ( k) ; 
end 


%calculate deviation from linear fit 
for i = 1:720 

RSS(i) = ( ( (y ( 1 , i ) -u ( 1 , i) ) A 2 + (y(2, i) -u (2, i) ) A 2 + (y(3,i)- 

u ( 3 , i ) ) A 2 + (y (4, i) — u (4, i) ) A 2 ... 

+ (y (5, i) -u (5, i) ) A 2 ) / 4 ) A .5; 

end 
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subplot (2,1,1) ; 
plot (theta, WSS) 

xlabel ( ' \ theta [degrees] ' , ' FontSize ' , 20 ) 
ylabel ( ' \tau_{ wall } [N/m A 2 ] ', 'FontSize', 20) 
axis ( [min (theta) , max ( theta) , 0 , max (WSS ) ] ) ; 
grid on 

subplot (2,1,2) ; 
plot (theta, RSS) 

xlabel ( ' \ theta [degrees] ' , ' FontSize ' , 20 ) 
ylabel ('RSS error' , 'FontSize', 20) 

axis ( [min ( theta) , max (theta) , 0 , max (RSS ) ] ) 
grid on 
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Appendix D — Supplemental DVD 

A supplemental DVD is available for this report. The disc contains the data, code, figures and files 

that were necessary for the generation of this work. Also included is an electronic copy of this report. 

Each file will be explained such that future use will be as simple and straightforward as possible. The 

supplemental DVD can be obtained by contacting the NASA Center for Aerospace Information, 

7121 Standard Drive, Hanover, Maryland 21076-1320, phone: 301-621-0390. 

D.l The “Unidirectional” Folder 

The contents of this folder include the following: 

1 . A tile titled “Summary RE l 7700 5 inch.” This is a summary of the unidirectional case in which S 
= 5 in. (127 mm), Re = 17700. Presented in this file are ensemble average velocities across the tube (a 
single profile) and across the disc spaces (five profiles). Also presented are 77 plots in the disc and 
tube spaces). 

2. A file titled “Summary_RE_7600_2. 125 inch.” This is a summary of the unidirectional case in which 
S= 2.125 in. (54 mm), Re = 7600. Presented in this file are ensemble average velocities across the 
tube (a single profile) and across the disc spaces (five profiles). Also presented are TI plots in the disc 
and tube spaces). 

3. A file titled “Summary_RE_7600_5_inch.” This is a summary of the unidirectional case in which S = 
5 in. (127 mm), Re = 7600. Presented in this file are ensemble average velocities across the tube (a 
single profile) and across the disc spaces (five profiles). Also presented are TI plots in the disc and 
tube spaces). 

4. A file titled “Summary RE l 7700 2. 125_inch.” This is a summaiy of the unidirectional case in 
which S = 2.125 in. (54 mm), Re = 17700. Presented in this file are ensemble average velocities 
across the tube (a single profile) and across the disc spaces (five profiles). Also presented are TI plots 
in the disc and tube spaces). 

5. A folder titled “Uni flow_viz_pics.” In this folder there are 4 files. These files are: A file titled 
“5high.jpg.” This is a flow visualization result for the case Re = 17700, S = 5 in. (127 mm). A file 
titled “five_low.jpg.” This is a flow visualization result for the case Re = 7600, S = 5 in. (127 mm). A 
file titled “twohigh.jpg.” This is a flow visualization result for the case Re = 17700, S = 2.125 in. 

(127 mm). A file titled “twolow.jpg.” This is a flow visualization result for the case Re = 7600, 

S = 2.125 in. (127 mm). 

6. A folder titled “RE_17700_2.125_inch.” In this folder, there are 10 files called “rl.jpg,” “r4.jpg,” 
“r5.jpg,” “r7.jpg,” “requalsl.xls,” “requals4.xls,” “requals5.xls,” “requals7.xls,” “requals8.xls” and 
“tube_RE_17700_2.125_inch.” The 5 jpg files are pictures of time and FFT traces for the radial 
locations r/D = 0.1 18, 0.473, 0.591, 0.827, and 0.945 (1, 4, 5, 7, and 8 in.). The Excel (Microsoft 
Corporation) file bearing the name “requalsl.xls” is time-history data sets for the radial location r/D = 
0.1 18 (r = 1 in.). Included in this Excel file is the raw voltage data that the anemometer recorded on 
the sheet “voltage, ” the converted velocity values and the ensemble averaged velocities and 77’ s on 
the sheet “velocity. ” Also provided in this file are plots of average velocity and TI. Finally there is a 
sheet called “spectral” in which there is a single, high frequency sample that was used for FFT 
purposes (for reference, this spectral sample can be processed with the code xlsfft.m, as described in 
appendix C section C.7). The other excel files are all similar in nature to the one that was just 
described. The exception is the file “tube_RE_17700_2.125_inch,” in which there is a single profile 
in the tube space. 

7. A folder titled “RE_1 7700_5_inch.” The description of the contents of this folder is identical to that 
of number “6” above. 

8. A folder titled “RE 7600 2. 125 inch.” The description of the contents of this folder is identical to 
that of number “6” above. 
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9. A folder titled “RE_7600_5_inch.” The description of the contents of this folder is identical to that of 
number “6” above. 

D.2 The “Oscillatory” Folder 

The contents of this folder include the following: 

1 . A folder called “flowmovies.” Within the flow movies folder there are three Microsoft Power Point 
files called: “30_RPM_2.125_spacing.ppt,” “30_RPM_5_spacing.ppt” and 

“70_RPM_2.125_spacing.ppt.” These are “movies” of oscillatory flow visualization results compared 
to computational work done at CSU. The titles of the files are explanatory of the contents. 

2. A folder called “2.125-30.” This is a set of data for the operating condition to = 30 rpm, S = 2.125 in. 
This folder contains seven sub-folders called “ssl,” “ss2,” “ss3,” “ss4,” “ss5,” “ss6,” and “tube.” All 
of these sub-folders contain similar information, so they will all be described at once. To begin, in 
each sub folder there will be an ASCII file that has to do with probe position. This file will be a 
column vector of probe sampling positions (more than ten locations and less than 100). All of the 
positions are in inches. Along with the position file there will be a collection of ASCII files with 
numerical titles. These titles correspond to the probe location at sampling time. Within these files are 
two column vectors, one of the columns is the velocity signal (as delivered from the anemometer) and 
the other column is a signal from the TDC indicator. This data is a continuous sample of the flow and 
the TDC indicator for 100 cycles. TDC is marked by a value near zero, away from TDC, the value is 
approximately 4. Also in each of these folders are two files called “pTI” and “pVel.” This stands for 
“processed” 77 or velocity. In other words, those files are the ensemble-averaged velocity and 7/ data 
for a entire cycle at each radial location. Using the code countV.c (as described in appendix C section 
3.6), the position file and the collection of velocity /TDC signal files, the files pTI and pVel can be 
generated. The general layout of the output files is as follows: 



position 1 

position 2 

position n 

o 

in 

o 

u(0.5°,positon 1) 

u(0.5°,positon 2) 

u(0.5°,positon n) 

i° 

u(l°,positon 1) 

u(l°,positon 2) 

«(1 ° ,positon n ) 

360° 

u(360°,positon 1) 

u(360° , posit on 2) 

u(360°,positon n) 


3. A folder called “2.125-70.” The contents of this folder are identical in nature to that described in 
number 2. The difference is that the operating frequency was 70 rpm and the disc spacing was 
2.125 in. 

4. A folder called “5-30.” The contents of this folder are identical in nature to that described in number 
2. The difference is that the operating frequency was 30 rpm and the disc spacing was 5 in. 

A couple of notes for the future user: The above data (e.g., pTI or pVel) can be plotted with the 
Matlab routine discthreed.m as described in appendix C section C.8. 

D.3 The “Codes” Folder 

In this folder are the following codes: 

1) AC code called “stepper_program.c” as described in appendix C section C.2 
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2) AC code called “dc.c” - as described in appendix C section C.3 

3) A Matlab code called “hotwirechambercalibration.m” as described in appendix C section C.4 

4) AC code called “oschwtraverse.c” as described in appendix C section C.5 

5) AC code called “countV.c” as described in appendix C section C.6 

6) A Matlab code called “xlsfft.m” as described in appendix C section C.7 

7) A Matlab code called “discthreed.m” as described in appendix C section C.8 

8) A Matlab code called “oscwallss.m” as described in appendix C section C.9 

There is also a file called “cib.o.” This file is a collection of commands necessaiy for using an IEEE 
488 bus/card with the above described programs. This file has to be included at the command line for 
proper compellation. For example, the correct syntax for running the code dc.c is as follows 

gcc -o ex cc u tab I ell I cn a me dc.c cib.o -lm 

D.4 The “Miscellaneous” Folder 

This is a collection of figures and so forth that went into creating this report. We will not try to 
describe them all. If the reader is inclined to satisfy his or her curiosity then so be it. 
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